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SUMMARY 


A procedure for calculating three-dimensional, compressible laminar boundary-layer 
flow on general fuselage shapes is described. The boundary-layer solutions can be obtained 
in either nonorthogonal body-oriented coordinates or orthogonal streamline coordinates. 
The numerical procedure is second-order accurate, efficient and independent of the cross- 
flow velocity direction. 

Numerical results are presented for several test cases, including a sharp cone, an el- 
lipsoid of revolution, and a general aircraft fuselage at angle of attack. Comparisons are 
made between numerical results obtained using nonorthogonal curvilinear body-oriented 
coordinates and streamline coordinates. A user’s manual with a detailed description of 
computer programs and input is presented in Volume II. 
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NOMENCLATURE 


A, B stagnation point velocity gradients 

ft major and minor semiaxis lengths of the ellipsoid of revolution 

a 0 cylinder radius 

C PnlptV-t 

C* BjA 

c p*!p 

G ' fx skin friction coefficient in the x-direction based on the edge condition, 

Eq. (103a) or Eq. (104a) 

G/y skin friction coefficient in the y-direction based on the edge condition, 

Eq. (103b) or Eq. (104b) 

Cp Pressure coefficient 

c p specific heat 

E H/H e , Eq. (48) 

F ufu t , Eq. (48) 

fi F, Eq. (48) 

G v/V ref or v y /V ref , Eq. (48) or Eq. (56) 

9( G, Eq. (48) or Eq. (56) 

H total enthalpy 

hi,h ,2 metric coefficients in the x and y coordinates, respectively. 

*\ 3i k indices in the x, y, and z direction, respectively 
i max, jmax , kmax 

number of boundary-layer grids in the x, y, and f direction, respectively 
K coefficient of thermal conductivity (= c p p/Pr) 

Ki,K 2 geodesic curvature of the curves y = const, and x = const., respectively, 
Eq. (5) or (19) 

Ku,K 2 i parameters defined in Eq. (6) or (19) 
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M 

ni,n 2 ,n 3 

P 

Pr 

r 

R, 0, <f> 
Rc oo 
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Mach number 

coefficients, defined in Eq. (54) or Eq. (60) 
coefficients defined in Eq. (54) or Eq. (60) 
pressure 

Prandtl number (0.72 ) 

radius measured from the X axes, Fig. 41 

spherical polar coordinates, Fig. 41 

free stream Reynolds Number, pooVooa/ Poo 

heat transfer at the wall, Eq. (108) 

arc length measured along y — const lines. 


T temperature 

t b/a 

u,v,w velocity components in the x,y, and z directions 

ujj,ue,u^ inviscid velocity components in the R,Q,<f> directions 

u x >,Uyi,u z i inviscid velocity components in the x', y 1 and z' directions 

velocity components in the x*,y *, and z* directions (near the stagnation point) 

dvjdy 

total velocity, Eq. (7) 

body-oriented coordinates (Fig. 1) or streamline coordinates (Fig. 2) 
rectangular coordinates with the origin at the nose point (Fig. 41) 
rectangular coordinates with the origin at the stagnation point, Fig. 37 or 38 
axial distance measured from the nose, see Fig. 1 
angle of attack 

Ax, Ay, Af grid spacing in the x, y, f directions, respectively. 

8 boundary-layer thickness; {z)v/v c = 0.995 

8* displacement thickness, defined in Eq. (107) 

e small angle to locate the initial streamlines near the stagnation point, Fig. 41 


u*,v*,w 


Vye 

V 

x,y,z 


x',y',z' 


x*,y*,z* 


X 

a 


via 


f transformed normal coordinate, Eq. (49) 

0 angle between x and y coordinates 

0 C half cone angle, Fig. 8 

0 r angle between two coordinate systems, ( x',y',z ') and (x*,y ,z), Fig. 37 

/j, molecular viscosity 

v n/p 

p density 

<f> azimuthal angle, 0 and it on the windward and leeward plane of symmetry, 

respectively, see Fig. 1 


subscript 


aw 

adiabatic wall 

b 

body-oriented coordinates 

€ 

edge of the boundary-layer 

osp 

origin of spherical polar coordinates 

S 

stagnation point 

st 

streamline coordinates 

t 

total 

w 

wall 

X 

partial differentiation with respect to x 

y 

partial differentiation with respect to y 

c 

partial differentiation with respect to £ 

oo 

free stream 
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1. INTRODUCTION 


Three-dimensional boundary-layer flows have been numerically studied for over three 
decades. During this period the capability to obtain numerical solutions has advanced 
from solving the similarity equations for relatively simple geometric shapes to the full 
nonsimilar equations for more complex configurations. The earliest referenceable numerical 
work, to the authors knowledge, was that of Raetz [1] and Der and Raetz [2]. These early 
papers remain as important contributions in that they introduced the stability of the 
mixed parabolic-hyperbolic system associated with the governing equations; i.e., the zone 
of influence-dependence principle. Blottner [3] presented a state-of-the art review of three- 
dimensional boundary-layer procedures that, with the exception of recently developed 
numerical methods, remains current at the present date. More recent treatments of the 
subject are presented in References [4] and [5]. 

Over the past decade, the major emphasis in computational fluid mechanics has focused 
on numerically solving the Euler and Navier-Stokes equations for increasingly more com- 
plex aerodynamic shapes. In many instances the Navier-Stokes approach is the only viable 
procedure, i.e., for flows with strong interaction and separation. However, Navier-Stokes 
solutions are generally much more expensive to obtain in terms of computer resources than 
boundary-layer procedures, and while capable of simulating the physics of complex flows 
they are often of low resolution due to grid point restrictions. Furthermore, Navier-Stokes 
solutions are not essential for many design and analysis procedures. 

Renewed emphasis on drag reduction [6], laminar flow control [7], and transition pre- 
diction [6] for complex flight configurations has clearly indicated the need to develop three- 
dimensional boundary-layer software that can be routinely applied to aerospace vehicles at 
a fraction of the cost associated with solutions obtained from the thin-layer Navier-Stokes 
equations. Research at the NASA Langley Research Center has resulted in the develop- 
ment and verification of two robust boundary-layer procedures for application to general 
aerospace configurations. One of these, a fourth-order accurate procedure for solving the 
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three-dimensional boundary-layer equations for aerospace configurations, was recently re- 
ported in Reference [8], 

Theory and equations for the other, a second-order accurate finite-difference procedure 
independent of the sign of crossflow velocity component for solving the three-dimensional, 
compressible, laminar boundary-layer equations, are presented in the present report. The 
software used to generate the numerical results has been optimized for fuselage shapes 
having a plane of symmetry; however, the general procedure is not limited to fuselage 
shapes and has been applied to wing flows [9], 

Results are presented for several test cases including a circular cone, an ellipsoid of 
revolution and a general aircraft fuselage at angle of attack. The method is valid for 
perfect gas flows from subsonic to hypersonic Mach numbers. Interaction between the 
inviscid and viscous flow is not included. A detailed description of the software including 
input/output for a fuselage shape is presented in volume II. 
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2. COORDINATE SYSTEM 


Non-orthogonal curvilinear surface coordinates are the most general system for the 
boundary-layer equations for complex aerospace configurations (see Ref. [5].). Two co- 
ordinate systems are presented in the present report: (l) a nonorthogonal body-oriented 
coordinate system with cross-flow planes perpendicular to the body axis (Fig. 1); and (2) 
an orthogonal streamline coordinate system (Fig. 2). 

Each of the two selected coordinate systems has its particular advantages and disad- 
vantages. The nonorthogonal body-oriented system is optimum from the viewpoints of 
grid generation and grid spacing control. Also, in certain aspects the interface software 
is simpler to apply since most inviscid solutions for bodies having a plane of symmetry 
use one coordinate plane perpendicular to the body axis. However, the boundary-layer 
equations are singular at the nose of the body (X = 0) and either a special transformation 
such as that used in Reference [10] or other procedures must be used to isolate this sin- 
gular point. The streamline coordinate system is orthogonal with zero values of cross-flow 
velocity at the wall and edge boundaries. The system’s origin is located at the stagnation 
point and is free of geometric singularities. But, the system is not independent of angle 
of attack, and downstream grid line distribution and grid point spacing is difficult, if not 
impossible, to control without an adaptive grid procedure such as that used in Ref. [11]. 
The primary interest in the streamline coordinate system in the present paper is in the 
eventual application of the software package to transition prediction; i.e., the output along 
the streamlines will serve as input for transition prediction procedures [12]. 
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3. THREE-DIMENSIONAL BOUNDARY-LAYER EQUATIONS 


3.1 Body-Oriented Coordinate System 


The governing equations (dimensional) for the body-oriented coordinate system are as 
follows: 

continuity equation 


c? d d 

— — [puh 2 sin0) + — (pvhi sind) + — — (pwhih 2 sin 0) = 0 
dx ay dz 


( 1 ) 


x-momentum equation 


pu du pv du du , __ , . „ 

— — — h — — — h pw — pu K x cot 0 + pv K 2 esc 0 + puvK 12 

hi dx hi dy oz 

esc 2 0 dp cot 0 esc 0 dp d . du . 


hi dx 


dy dz dz J 


( 2 ) 


y-momentum equation 


pu dv pv dv dv 2 r*- . 2r . * 

— — — — h pw— pv K 2 cot 0 + pu K\ esc 0 + puvK 2i 

hi dx hi dy dz 

cot 0 esc 0 dp esc 2 & dp ( d f dv ^ 

hi dx h 2 dy dz ^ dz 

energy equation 

pu dH pv dH 


( 3 ) 


dH 

hi dx h 2 dy ^ dz 


.L[. 

dz \ 


p dH 


+ M 1 


Pr’dz K 2 J j 


( 4 ) 


Pr dz 

The metric coefficients hi and h 2 are functions of x and y. The parameters K\ and K 2 are 
geodesic curvatures of the curves y = const and x = const respectively, where 

dh 2 


K ‘ = {^<^ cos9 ) - ^r} > K ' = - d ~t) 


( 5 ) 


and 


K n = 


hih 2 sin 0 


21 hih 2 sin 2 0 { 


/ o d hv i 

(1+cos 


dh 2 

~d^ 


(1 + cos 2 0) 


2cos<? 


2 cos 0 


dx 

dh 


) 

?} 


( 6 ) 
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V is the total velocity and is given by 


V = (u 2 + v 2 + 2 uv cos 9) 1/2 


( 7 ) 


The boundary conditions are 


z = 6, u = u e (x,y), v = v e (x,y), H = H e 


dH< 


0, u = v = 0, u; = w w , H = H w or = 0 


(8.a) 

(8.b) 


At the edge of the boundary-layer, the pressure gradients are related to the inviscid veloc- 
ities by the following equations: 


fu* 

l hi 


Pe , + h . _ U \K X cot 6 + v]K 2 esc 6 + U'V'Ku ) 

ax h 2 dy ) 


esc 9 dp cot 9 esc 9 dp 
hi dx h 2 dy 


(9-a) 


Pe 


! + cot 9 + u]Ki esc 9 + u e v e K 21 ) 

[ hi dx h 2 dy ) 


cot 9 esc 9 dp esc 2 9 dp 


(9.b) 


hi dx h 2 dy 

The perfect gas equation of state and Sutherland’s viscosity are used to close the equation 
set. 


For the windward and leeward plane of symmetry t; and Ki esc 9 are zero and 6 is 
generally 7r/2 ( 9 is retained for the general system.). Consequently, each term in the y- 
momentum equation vanishes. However, partial differentiation of Eq. (3) with respect to 
y yields an equation for dv/dy. After differentiation and using the appropriate symmetry 
conditions ( du/dy = dw/dy = d 2 v/dy 2 = dH/dy = dhijdy = dh 2 /dy = 0) along with 
Eq. (9), the governing equations for the plane of symmetry become 
continuity equation 


Q 

—(puh 2 sin 9) + pv v hi 
dx 


(10) 


5 



x-momentum equation 


pudu du o „ .u.du, , _ r d . du. 

T l ai + l,w a-z- f ’ u K ' cote = ~ u '*‘ coU ’ ) + a^-g;) 


y-momentum equation 
pu dv y 


dv 


Pe( 


+ pw + -f- vl + puv v K 2 1 + pu 2 

~ fl2 

V + ■J~ + UeWye-f^2l) + 


/ll dx d,2 /i-2 

/ii dx h 2 


d(Ki esc 0) 

dy 


d(Ki esc 0) d . dv y . 
dy ^ dz ^ dz 


energy equation 


pu d# 


+ piw- 


dtf 


= -( 

dz \ 


H dH 

P^~dz +fi ^ Pr'dz K 2 


hi dx dz 

where v y = dvjdy, v ye = dv e jdy and V — u along these lines 
The boundary conditions for the plane of symmetry are 


1 )££>} 


in) 


( 12 ) 


(13) 


z = S, 


u 


= u e (x,y), v v = v ye , H = H e 


,dH s 


(14.a) 


z — 0, u — v = Vy — 0, w = w w , H — H w or (— — )„, = 0 (14. b) 

dy 


3.2 Streamline Coordinate System 


The streamline coordinate system is orthogonal; consequently, with the exception of 
certain metric coefficients and the boundary conditions, the governing equations can be 
obtained directly by equating 6 — 7t/2 in Eqs. (l)-(4) and Eqs. (10)-(13), i.e., 
continuity equation 


d d d 

— (puh 2 ) + —(pvhi) + — (pwh^i) = 0 


(15) 


x-momentum equation 


pu du 
h\ dx 


pv du du , _ _ 1 dp 

+ — — + pw— + pv 2 K 2 + puvK n = - — — 
/i 2 dy dz h\ dx 


d , du. 


( 16 ) 
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y-momentum equation 


1 dp d . dv. 


puSv pvdv dv , =~TJT + 

h ! dx h 2 dy F dz K h 2 dy dz dz 


(17) 


energy equation 

pu dH pv dH dH _ d j p, dH 
7 H~dx + T 2 ~dy +f>W dz ~ dz { Pr dz 

The parameters K u K 2 , K n , and K 2i are given by 

1 dh x „ 1 dh 2 


+ i ^ I" ~ "aT \ SZ ^ p r * dz' 2 J 


(18) 


K x = - 


K 12 = 


hih 2 dy 
1 dhi 


K 2 


h\h 2 dx 


hih 2 dy 

V is the total velocity and is given by 

V = (u 2 + v 2 ) 1/2 


1 dh 2 

= -K u K 2 \ = — — = K 2 


h\h 2 dx 


(I9.a) 

(I9.b) 


( 20 ) 


.dH , 

u = v = 0, w = w w , H = H w or (— J™ = 0 


(21. a) 
(21.b) 


The boundary conditions are 

z ■ — S) ll — Ug (j, y) j v 0, H H t 

z — 0, 

At the edge of the boundary-layer in this coordinate system, the pressure gradients are 
related to the inviscid velocities by the following equations: 

(22. a) 
(22. b) 


Pt 


u e du e 


1 dp 

hi dx hi dx 
p e u 2 e dhi 1 dp 


h y h 2 dy h 2 dy 

In the streamline coordinate system for the boundary-layer, the metric coefficient hi is 
defined as 


hi = 


(23) 


u. 


The governing equations for the plane of symmetry become 
continuity equation 

d ^ 

— (puh 2 ) + pVyhi + —{pwhih 2 ) = 0 
dx oz 


(24) 
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x-momentum equation 


pu du du u e du e d . du. 

+ pw ~d~ z = pt hi~fa + -d~z^irJ 


y - momentum equation 


pu dv v 

x 

energy equation 


dv u 


*T57 + ^ aT + C' + fu ' , ' K " + ^ — a y • a, - a. 


-diiTi -d/fi 5 , dv„. 
2 — -p e „2_i + ( M ^i) 


pu dH dH 

+ pu) 


= JLiJt 

dz \P» 


hi dx dz 

where v y — dvjdy and V = u along these lines. 

The boundary conditions for the plane of symmetry are 


dH „ 1 . d ,V\ 

Pr dz + ^ Pr^ dz 2 


(25) 


(26) 


(27) 


2 = (5, u-u e (x,y), Vy = 0, H — H t 
z = 0 . 


,dH , 


u = v = u„ = 0, w = w w , H = H w or (-r— = 0 

dy 


(28. a) 
(28. b) 


3.3 Three-Dimensional Stagnation Point 


To obtain the boundary-layer solutions at the three-dimensional stagnation point, the 
governing equations for three-dimensional laminar compressible flows in Cartesian coor- 
dinates are required and can be obtained by setting hi = 1, h 2 = 1> and 0 = 7r/2 in 
Eqs. (l)-(4). Superscript * is used to distinguish this coordinate system from the other 
coordinate system, 
continuity equation 


d 

dx* 


(pu*) + 


d 

dy* 


(pv*) + 


d 

dz* 


(pw*) = 0 


(29) 


x-momentum equation 


(30) 
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y-momentum equation 


t dv* t dv* dv* dp d , dv * . 

a^ + l>v W + pw a7 = ~S7 + ^aP ) 


(31) 


energy equation 


t dH dH t dH d 

pu 'a ? ' + pv ai- + pw —■ ~ 


dz* dz* Pr dz 


\- 

\ Pr 


dH . 1 , d .V* 2 ' 


(32) 


where 


V* = (u* 2 + v* 2 ) 1/2 


(33) 


The boundary conditions are 


z* = 8, u* = u*(x,y), v* = v* t {x*,y*), H = H e 


z* = 0, u* = v* = 0, 


3 ff 

H = H w or (— ) w =0 (34. b) 


(34. a) 


3.4 Sharp Cone 


The boundary-layer solutions for the flow on the sharp cone are used for generating 
initial profiles near the nose tip for sharp-nose fuselage shapes. The governing equations 
for a sharp cone, written in terms of x, the coordinate along a cone generator, y, the 
cone azimuthal angle and z, the coordinate normal to the cone surface can be obtained by 
substituting h x = 1, h 2 = r = xsin0 c , and 0 = tt/2 from Eqs. (l)-(4) and Eq. (10)-(13) 
where 0 e is the half cone angle: 
continuity equation 


3 3 3 

—(pux sin0 c ) + — (pi;) + — (ptuxsin 0 C ) = 0 


(35) 


x-momentum equation 


du pv du du p , d , du. 

pu w ~ + — nrr + pw * v = ) 

dx xsin 6 c dy dz x dz dz 


(36) 
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y-momentum equation 


dv i ov dv dv puv 

pu 1 — — — h pu >— — I 

dx zsin 0 C dy dz x 


1 dp 
x sin 6 C dy 


d dv 

+ Tz^ai ] 


(37) 


energy equation 
dH 


pu 


dx 


+ 


dH dH 

b pu>— — 

xsin0 c dy dz 


pv 




where 

V = (u J + 


(39) 


The boundary conditions are 


z = 6, u = u e (y), v = v e {y), H = H e (40. a) 

dH 

2 = 0, u = v = 0, w = w w , H = II w or (— )«, = 0 (40. b) 

The conical inviscid flow assumption has been made in the above equations, i.e., all gradi- 
ents of the inviscid variables in the x direction are assumed to be zero. 

At the edge of the boundary-layer, the pressure gradient in the y-direction ( dpjdy ) is 
related to the inviscid velocities by the following equation: 


p e v e dv e , 1 dp 

+ p e u e v e = — : 


(41) 


sin0 c dy ' sin 0 dy 

Using the conical inviscid flow assumption, the following equation can be obtained from 
Eq. (9. a). 


1 du. 


= v e 


sin 6 C dy 

The governing equations for the plane of symmetry become 
continuity equation 

— (pux sin0 c ) + pv v + —(pwxsmOc) = 0 
dx oz 


(42) 


(43) 


z-momentum equation 


du du d . du. 

+ pw Tz = Tz^ 


( 44 ) 
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y-momentum equation 


dv v 

pu— h pw 

ox 


dvy 

dz 


+ 


+ puVy 

x sin 9c dy x 


= Pc{ 


V l* , , d ( „ dv V} 

xsin#,. x dz dz 


(45) 


energy equation 


dH 

pu— + pw 

ox 


dH 

dz 


d_ 

dz 


H dH 
Pr dz 


+ ^(1 


Pr’dz y 2 


>) 


(46) 


where v v = dv/dy, and V = u along these lines. 

The boundary conditions for the plane of symmetry are 



u = u e (y), v y = v ye , H = H e 
u = v = v u = 0, w = w w , H — H w or 


,dH . 

( ^ = 0 


(47.a) 

(47.b) 
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4. TRANSFORMED EQUATIONS 


The boundary-layer equations are transformed to the coordinate system used by Cebeci 
et al. [13] which removes the singularity at x = 0 and allows coupled solution of the 
continuity and momentum equations. 

The following definitions are introduced 


F = ft = uju e , G = g { = v/V ref , E = H/H e 


(48) 


together with the transformation 

x = x, y=y, 

where 


= ML r t dz 

V fi e s Jo p e 


(49) 


s = J hidx (50) 

In the present report, V re f is chosen to be Vqo except for the three-dimensional stagnation 
point equations, where it is taken to be v e . 


4.1 Body-oriented coordinate system 

Using the transformation given in Eq. (49) and the relations given by Eq. (9), the 


governing Eqs. (l)-(4) are transformed to the following form: 
x-momentum equation 

F = ft (51. a) 

(CF { ) S + mifF ( — m-iF 2 — m^FG + m 6 F ( g — m 8 G 2 + mnc — m 13 .F ? 

= mio(FF x — F ( f x ) + mj (GF V — F ( g v ) (51. b) 

y-momentum equation 

G = g ( (52. a) 

(C^Cr^); -)- mifG ( — m 8 G 2 — m^FG + m^G^g — trigF 2 -|- mi2C — ixiaGf 

— m io {FG X — G ( f x ) + m. 7 (GG y — G ( g u ) (52. b) 
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energy equation 


[niE ( ) ( + n 2 E { + (n 3 ) f - m 13 £ ( = m 10 (FE x - EJ X ) + m 7 (GE v - E ( g v ) (53) 


The coefficients mi to m 13 and n \ to n 3 are as follows: 




s du e ) s d 


mi hiu e dx \ hih 2 sin 6y/pJJT e dx 

s du e 


{h 2 sin Oy/pT^l} 


m 2 


— sK\ cot 6 


h\u e dx 

V„, , av„, 

u e h 2 u e ay 

s dV ref 


m 4 = sK 2 \ + 


hiV ref dx 


m 5 


m 6 


m 7 


sV ref du e V ref 

T r 

h 2 u l e dy u. 


hih 2 s\n6 s fp^r e 

sVrcf 

h 2 u e 

m 3 = sK 2 cscO (—^) 2 
u. 




mg = sK\ esc 6 


m w = — 
« 1 


u. 


f rtJ 


I 1 du e u e du e , 

mu = 5 t - — + - — 7— cot dKi + CSC 0 K 2 { 


m 12 


\hiu e dx h 2 ul dy 
s [ u e dv e ^ v e dv t 


m 13 


n 1 


u e V ref dx h 2 dy 
_ (pw) w [P^S 


U e U e ) 

— cot 6K 2 v l + esc 9K\u\ + K 2 iu e v e j 


p e u t 


Me 


Pr 

n 2 = mif + m & g 

ns = ^(1 - i) {FF S + GG^) 1 + ^j-co S e(FG, + f,G)} 

The boundary conditions are 

f = 0: f = F = g = G = 0, E' = 0 or E = E w 
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(54. a) 
(54. b) 
(54.c) 
(54. d) 
(54.e) 
(54. f) 
(54.g) 

(54. h) 

(54. 1) 
(54.j) 
(54.k) 

(54.1) 

(54. m) 

(54. n) 
(54.o) 
(54. p) 

(55. a) 



£ = 6 : 


F= 1, G = v e /V ref , E = 1 


(55. b) 


The governing equations for the plane of symmetry are transformed by defining 


F — = u/u e , G = g { = v v jV rc f, E = H/H e (56) 

x-momentum equation 

F = f t (57. a) 

(CF f ) f + mifF { - m 2 F 2 + m 6 F ( g + m n c - m i3 F ( = m l0 (FF x - F ( f x ) (57.b) 

y-momentum equation 

G = j, (58.a) 

(CG t ), ■f wii/Gf — m 3 G 2 — 1714FG -f m 3 G^g — mgF 2 + mi 2 c — 

= m 10 {FG x - GJ X ) (58.b) 

energy equation 

{ n iF()s + n 2 F { + («3)f — mi 3 E ( = m,io(FE x — E ( f x ) (59) 


The coefficients mi,m 2 ,m 4 ,mi 0 ,m 1 3 ,ni,n 2 are the same as in Eq. (54). The remaining 
coefficients are defined as follows: 


sV, 


m 3 = 


ref 


h 2 u e 

m e = m 3 

su e d{K\ esc 6 ) 


m 9 


mu 


m 12 


Kef dy 
s du. 


hiu e dx 
su e d{K\ esc 6 ) 


s cot 6 K 1 


Vref 

Cu 2 


dy 


s , 1 dv, 

+ — ( 


ye 


V rc / hi dx U t h<i 


+ Hr + K n v ve ) 


”3 = ^(1 - -bz)FF { 


Pr 


(60. a) 
(60. b) 
(60.c) 

(60. d) 

(60.e) 

(60.f) 
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1 I S 

mi = 2 l 1 + S7u, 

5 du e 
m2 h\u t dx 

m 3 = 0 

m 4 = 5-^21 = —sK 2 
m 5 = 0 
m 6 = 

5Voo 


- — } + 3 =JT {^a y/pTi^l} 

l e dx J hihiy/pTJ^e dx 


The boundary conditions are 

£ _ q : f = F = g = G = 0, 22' = 0 or E = E w (61.a) 

f-fc! r=l, G = ^, i = l (61 - b) 

4 .2 Streamline Coordinate System 

The transformed equations have the same form as Eq. (51)-(53) off the plane of sym- 
metry. The coeffcients m, to m„ and n, to n„ obtained by setting « = */2, V„, = V„, 
^ = uJVoo, and v e = 0 from the Eq. (54), are as follows: 

" (62. a) 

(62.b) 
(62. c) 
(62. d) 
(62 .e) 
(62. f) 

(62. g) 

(62.h) 

(62. 1) 
(62. j) 

(62.k) 

(62.1) 

(62. m) 


u 

r , 

hih 2 yfpjl^u^s dy l «e J 


r«7 = 


h 2 u e 


mg = sX 2 (— ) 2 
u c 

mg = sKi — 

V oo 
5 

mio = — 

/ii 


m u 

s du e 

H 

JS 

1 

mi2 

u, 

= sK 'y- 

v oo 


{pw)v> [ 

mi3 



p e u e V 




ni = p^ 

n 2 = mi/ + m 6 ff 

Cu 2 


n 3 =^(f + 
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(62.n) 
(62. o) 
(62.p) 



The boundary conditions are 


f = 0: f = F = g = G = 0, E' = 0 or E = E w (63.a) 

f = & : ^=1. G = 0, E = 1 (63. b) 

The transformed equations for the plane of symmetry for this coordinate system has 
the same form as Eq. (57)-(59). The coefficients m 1 ,m 2 ,m 4 ,m 1 o,m 13 ,n 1 ,n 2 are the same 
as in Eq. (62) . The remaining coefficients are defined as follows: 


tv,, 

m 3 = 

h 2 u e 

(64. a) 

m 6 = m 3 

(64. b) 

su e dKi 

m g — 

Voo dy 

(64. c) 

s du, 

m,, - 

h\U e dx 

(64. d) 

su e dKi 

mi 2 = 

Vco dy 

(64. e) 

n * = C >-^ FF < 

(64. f) 


The boundary conditions are 


f = 0: f = F = g = G = 0, E' = 0 or E = E w (65.a) 

( = (•• F= 1, G — 0, E = 1 (65. b) 
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4.3 Three-Dimensional Stagnation Point 


The governing equations for three-dimensional laminar compressible flows in rectangu- 
lar coordinates, Eq. (29)-(32), are transformed using Eq. (49) and Howarth’s [14] inviscid 
velocity components near the stagnation point. For s approaching zero, the following sim- 
ilarity (ordinary differential) equations are obtained, 
x-momentum equation 

(C/7 + //" - (/')’ + | f"g + ^ = o (66) 

y-momentum equation 

(C 9 ")' + fg" - |(s') ! + jg"g + jj=° (67) 

energy equation 

(§;&)' + (/ + 1 g)E' = o (68) 

The equations above are based on the assumptions that the outer flow is irrotational and 
that the inviscid velocity components near the stagnation point can be approximated by 


u. 


Ax*, 


v* e = By* 


(69) 


Equations (66)-(68) can be obtained from Eqs. (5l)-(54) by substituting u* t = Ax* , V Tt j = 
v; = By*, h . l = i,h 2 = i,9 = tt/2, and taking limit as s approaching to zero. The primes 
denote ordinary differentiation with respect to f, i.e., 


f , = dj_ 
1 d$ 


u 

u! 


, _ dg 
9 d f 


v* , „ H 

— , and E = —— 

K H e 


(70) 


The boundary conditions are 

g — 0 : f = f' — 9 = g' = 0, E 1 — 0 or E = E w (71. a) 

f = /' = !, 9' = h E = 1 (71 -b) 
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4.4 Sharp Cone 


Using Eq. (49), the boundary-layer equations off the lines of symmetry are transformed 


to (Here, V ref = V*, ): 
x-momentum equation 


(c/„) t + j//„ - 


r | / ^co g^(PeA^e) 1 VgVoo . 2 

— — /ffl'f + L . ~ a — =- )fi { g + ( — ) g ( 

U‘ 2sm 0 c u e p'He 2 u] u e f 


[pw)w l p<U'S . _ Vqo , < 9/ f dg 

p e u e V n e {{ sin 0 c u e 9t dy i{ dy 


(72) 


y-momentum equation 


, . 3 1 dv e 

( 0ff)f + - sin 0 c u e dy 


*? " 


(P e Me) 1 W.V'oo 


2 sin 0 c u e p e *t e 2 u\ 


)9 tS 9 


+(in 


dv e v e p e ( pw) w fp^s 

+ 77-)- ~ -rz~\—9« 


sintfetieVoo dy ’p p e u e 

energy quation 




Voo . dg ; dg 
sin 0 e u} 9i dy 9{( dy^ ^ ^ 


(§- r E i'h +11 / + ( „Z“.. I E < - — . ^E- 


+ 


P e Uc V 

dg. 


2 J ' '2 sin 6 c u e p c p, e 2 u* 

{|d - - ** 

where f { = df /d$ = u/u e , g ( = dg/dg = v/V^, and E = H/H e 
The boundary conditions are 


) 


(74) 


<r = 0: 

f = fc : 


/ = ft = 9 = 9( = 0 , E C = 0 or E = E v 


ft = 1. fc = ^-, ^ = 1 

^ oo 


The boundary-layer equations on the lines of symmetry become 
x-momentum equation 

(<?/«)* + + sinTu f * <9 = ° 

l sin t/ c u e 

y-momentum equation 

( C 9 S ()( + \f9it ~^T9 2 ( ~ f{9 ( + ■ V T 9 k9 + ( v . t>y *. v , 

° u.sm^ f sinS c u t t^F^sin#,. p 


. P* 

' t r * 


(75. a) 
(75. b) 


(76) 


= 0 (77) 


18 



energy equation 


+ + &<>} E ‘ + {w^~ h )C! ^<) ! - 0 < 78 ) 


where f ( = df/dg = u/u e , g t = dgjdg = vjV w and E = H/H t 
The boundary conditions are 


f = 0 : f = ft =9 = 9 t = 0, E f = 0 or £ = £„ (79.a) 

S = 6 : ft = 1» E = 1 (79. b) 

^OO 

Equations (72)-(74), and (76)-(78) can be obtained from the Eqs. (51)-(54) and Eqs. (57)- 
(60) by substuting hi = l y s = x, hj = xsin0 c , 6 = 7r/2, V re ^ = Vqo and with the conical 
inviscid flow assumption, i.e., dujdx = 0, dp/dx = 0, dpjdx = 0, dpjdx = 0. 
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5. NUMERICAL METHOD 
5.1 Blottner’s Iterative Method 


The transformed stagnation point equation, Eq. (66)-(68), the governing equations off 
the lines of symmetry on the sharp cone, Eq. (72)-(74), and the equations along the lines 
of symmetry on the cone, Eq. (76)-(78), are solved using Blottner’s iterative method [15]. 
All the equations listed above can be expressed in the following form: 
x-momentum equation 

F = f ( (80. a) 

(C'FfJf + mifF ) - m 2 F 2 - m 5 FG + m 6 F ( g - m s G 2 + m n c - m 13 F ( = m 7 {GF v - F ( g v ) 

(80. b) 


y- momentum equation 

G = g { 

(CG { ) { + mJG { - m 3 G 2 - m 4 FG + m 6 G { g - m 9 F 2 + m n c - m 13 G f = m 7 {GG y 


(81. a) 

G ( g v ) 

(81.b) 


energy equation 

(riiE ( ) ( + n 2 E ( + (n 3 ) f — mi 3 E ( = m 7 (GE y — E { g y ) (82) 

The above equations are linearized using Newton-Rhapson’s linearization technique [3]. 
The f-derivative terms are discretized using a central finite-difference scheme. To solve the 
governing equations off the line of symmetry of the cone (Eq. (72)-(74)), an implicit 
marching procedure ( Ref. [15]) is used. Here, for the y-derivative term, an implicit second 
order backward finite-difference is used. 

For abbreviation, finite-difference operators are defined as 

$ Ft = ^* +1 — ^ k ~- k = 2, 3, .., kmax — 1 (83. a) 

{ Af* + Affc-i 
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A f (C*A f F*) = 


[ r F k+1 -F k F k -F> 

\ C k+ 1/2 ZZ 1/2 

k-1 l 


Aft 


HhI 

Affc-i j 


W = 


Ac* + Ac* 

if j=2 

{A y 2 _ x - (Ayy-i + Ay ; - 2 ) 2 } Jy + (Ay,_i + Ayy_ 2 ) 2 ^y-i ~ Ay 2 _iEy- 2 


c r Ft- ^-1 


(83. b) 
(83.c) 


(At/j_i) 2 (Ay J _i + Ay ; _ 2 ) — Ayj_i(At/j_i + Ayy_ 2 ) 5 


if j > 2 


(83. d) 


where C fc±1/2 = \{C k + C k±1 ), Ac* = C*+i - ft. Ay, = yy+ 1 - Vy, and j and fc represent the 
y and C directions, respectively. The overlined quantity is the converged solution at the 
previous step (j — 1 or j — 2) . 

The finite-difference equations for the Eq. (80)-(82) can be written as follows: 
x-momentum equation 

ft ~ fk-i ~ —(ft + F „- .) = 0 < 84 ' a ) 

A,(C'i A t F t ) + + 6 s Fth - 7 l «;7*) - rn,(2FF t - 7 \) - m s (G t F t 

+F k G k - F k G k ) + m 6 (ff*<5 f F* + 8 ( F k g k - g k 8 { F k ) - m 8 (2GG* - G*) + m n c k - m 13 F s 
= m 7 (G k 8 y F k + 8 y F k G k — G k 8 y F k — 8 l F k 8 y g k — 8 i F k 8 y g k + 8 ( F k 8 y g k ) (84. b) 


y-momentum equation 

a - + c M ) = 0 («•») 

A { {C k A { G k ) + m x (J k 6 t G k + 8 ( G k f k - J k 8 ( G k ) - m 3 (2GG* - g\) - m A {G k F k 
+F k G k - F k G k ) + m 6 {g k 6;G k + 8 s G k g k - g k 8 ( G k ) - m g (2FF k - F\) + m n c k - m n G s 
= rm{G k 6 y G k + S y G k G k — G k 8 y G k — 8^G k 8 y g k — 8 i G k 8 y g k + 8 s G k 8 y g k ) (85. b) 


energy equation 

A f {n l>k A ( E k ) + n 2 k 8 ( E k + 8 t n a>k - m n E ; = m 7 {G k 8 v Ej - 8 { E k 8 y gj) (86) 

where the overlined quantities are evaluated from the previous iteration. The energy 

equation does not require linearization since it is solved after the momentum equations. 

21 


The finite-difference momentum equations, Eq. (84) and Eq. (85), are rearranged in 
2x2 block tridiagonal form as 


hk - hk ~ i •+ + -fffc-i) (87. a) 

-A k H k _ 1 + B k H k - C k H k+1 + a k h k = D k (87.b) 


where 



A k i B k , C ki a k are 2x2 matrices, and D k is a vector. 

These equations are solved by the Davis Modified Tridiagonal Algorithm (See Appendix 
A). The finite-difference energy equation, Eq. (86), is arranged into the linear tridiagonal 
matrix equation form as 


B kE k - 1 + D k E k -(- A k E k+1 = C k (88) 

where A k , B k , C k , and D k are scalars. 

This equation is solved using the Thomas Algorithm. The momentum equations and the 
energy equation are solved iteratively in a uncoupled manner until the converged solution 
is obtained. The converged solution is usually obtained within five iterations. 


5.2 Matsuno’s Finite Difference Method 

5.2.1 Formulation of Finite Difference Equations 

Matsuno s finite-difference method [16] is used to march the solution downstream 
from the initial data plane (velocity and temperature profiles specified at initial data 
plane; see Appendices B and C for detail). The method is a modification of the predictor- 
corrector form of the Crank-Nicolson scheme, which was originally suggested by Douglas 
and Jones [17] to apply to the three-dimensional boundary-layer problem. This scheme 
is half implicit in the f direction, explicit in the y-direction, noniterative and has second 
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order of accuracy [16]. Also, the scheme is highly vectorizable for computation because 
the crosswise derivatives are formed independent of the the sign of the crossflow velocity 
component. 

Again, for abbreviation, finite-difference operators are now defined as 


l 


_ 


i.fc-i 


Af* + Af*_i 




A ft + Aft 






F l« 


Aft 


-Cl 






Aft 


(89. a) 

2^|89.b) 


&v F U = 


(A yy-on^i.t - + (Ayy)»(f},» - *?-!.») 

AyyAy y _j(Ayy + Ayy_ a ) 


(89.c) 


where Cj t; tl / 2 — |(Cy >Jk + Cy,t±i)> an( ^ *>J» and & represents the x, y, and 2 directions, 
respectively. 

Figure 3 shows the finite-difference molecule for the scheme. To formulate the finite- 
difference equations for the predictor step, the nonderivative terms are given as the value 
at the previous step (i-th step), the x-derivative is obtained by backward differencing, 
the first derivatives of y and f are obtained by using central differencing at the previous 
step explicitly, and the second derivative of f is obtained using central differencing at the 
predictor step implicitly. For the corrector step, the nonderivative terms are given as the 
predictor values, the x-derivative is obtained by backward differencing, the first derivatives 
of y and f are obtained using central differences at the predictor step (t + 1/2) explicitly, 
and the second derivative of f is obtained by averaging the t'-th step and (* + l)-th step. 
The finite-difference equations which approximate Eq. (51) through (53) are formulated 
as follows: 

Predictor 

x-momentum equation 


r-i+l/2 ri+1/2 Affc— X f j^t+1/2 . p,t+l/2\ « 

Jj,k - /y,k-i 2 + ' — u 


(90. a) 
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A,(Ck A t j£ ,/2 ) + m, /;.* ^ - m 2 (^ k ) ! - m 6 fj , G*,* 

+ m e 4» ^^y.fc - + m„ 4» - m is' 5 <4.t 


£it+l/2 pt rt-M/Z ri 

= mio( ^‘ f, ‘ AX../2 ~ " Ax,/2 } + mr(G '-'‘ ~ *<*?,* ( 90b ) 


m+1/2 

* J i,k 


y-momentum equation 


4 1/! - *825 - ^(^ ,/! + iS!) = » 


(91. a) 


A({C' j k . A { G'^ 2 ) + mi f) k 6 { G' jk - m 3 (Gy t ) 2 - m 4 F' j k G' j k 
+m 6 g) k 6 { G* j>k - m g (F} k ) 2 + m 12 c) k - m n 6 ( G) k 


^i+1/2 

G J.fc 


P / t + 1 / 2 _ ft 

M j ^ ^ ^ 


f IP* J.* 

“ lo(f >'‘ Aii/2 ‘ Ax,/2 


-^G/.A4*) ( 9I - b ) 


energy equation 


Af ( n l;\jb ^ ) + n 2jf,fc + ^f n 3j,Jb m 13^^y,*: 


4i ,/! - % 




~ S,E).„ 


ri+1/2 

Jj,k 


fU 


Axi/2 


) + m 7 (G) ik 6 v E) tk - 6 ( E) >k 6 v g) ik ) (92) 


Corrector 

z-momentum equation 


M - " +1 - + 4t->) = 0 


fi+ 1 f»+l 
Jj.k ~ Jj,k - 1 “ 


(93. a) 


A t {g; 1 1/! A, ± ^ | + m, /£ l/! - m s f£ 1/j Gfc 1 ' 2 

tm .^ 2 i t F‘r /! - m 8 (G "‘' 2 ) 2 + m„c ‘«' 2 - m^jft 1 ' 2 = 


A Xi 


_ et ^fU_Ik ) + m7 ( G ;«/ 2 ^ ^ s ‘+ , /») 

y-momentum equation 


(93. b) 


- ?;V-i - + Gj.y_,) = 0 


(94. a) 
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,<+1/2 A {Gj+k + G ),k) 


l, A t 


+ m, S\T 6,G%'I 2 - m,(GZ l, r - m, F%‘“ GJJ 


,»+1/2a2 


-,.+ 1/2 ^,<+1/2 


{ Of * + l /^t 

+m 6 g%'!' S t G^‘ - ^(F^)' + m„ c£' ! - m^G^’ = 


Axj 


_ 6!C fgi* !)*_£*) + mi (G‘+ ,/! «„g£ 1/2 - « t G;y /J 


(94.b) 


energy equation 


v J„* +1 / 2 

1 n ij,* 




4+1 + E 


itij 


+ ««.f + { <<,r - m iA 


;+i/2 


,<+ 1/2 




3.* 


Ax,- 




1/2 /j'.fc 1 /;*,) 


Alj 


)+m 7 (G;^ ! «„4; ,/! 




,<+1/2 


W* 1/2 )( 95 ) 


where the superscripts t, t + l/2, and t'+l denote the i-th step, predictor step, and corrector 
step, respectively. 

Both the predictor and corrector finite-difference momentum equations (Eqs. 90, 91, 
93, 94) are rearranged in the 2x2 block tridiagonal form as Eq. (87) and solved by Davis 
Modified Tridiagonal Algorithm (See Appendix A). Each (predictor and corrector) finite- 
difference energy equation (Eqs. 92, 95) is arranged into the same linear tridiagonal matrix 
equation form as Eq. (88) and solved using the Thomas algorithm. Although there is cou- 
pling between the momentum equations and the energy equation, these equations can be 
solved in an uncoupled manner due to the quasi-linearization involved in the predictor and 
corrector scheme. 


5.2.2 Stability 

The mathematical character of the three-dimensional boundary-layer equations was 
shown by Raetz [l] to be hyperbolic in the x — y plane, resulting in the formulation of the 
zone of influence and dependence principle. The influence of the solution at any point is 
transferred by diffusion to all points on the line normal to the surface and by convection 
downstream along the streamline through that point. The zone of dependence for a certain 
point is a wedge shaped region facing upstream bounded by two characteristic surfaces each 
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containing the outermost streamlines ( one is the inviscid streamline and the other is the 
limiting streamline ) passing through the point. 

The zone of dependence designates the minimum amount of initial data to be sup- 
plied; in other words, the difference molecule must include the information in the zone of 
dependence. Because of this principle, the standard finite-difference methods for solving 
the three-dimensional boundary-layer equations using the numerical marching procedure 
in the x— and y-direction must be modified whenever the sign of the cross-flow reverses. 
More exactly, when the direction of any streamline in the boundary-layer is opposite to 
the numerical marching direction, a modified method must be used. The finite-difference 
methods used by Shevelev [18], Dwyer and Sanders [19], Mclean [20], and the Box scheme 
used by Cebeci et al. [13] axe examples of methods which require modification for the region 
where the crossflow direction is opposite to the numerical marching direction. The Zig-zag 
scheme used by Krause [21], Zig-zag Box scheme, and Characteristic Box scheme [10] are 
examples of modifications used in this region for standard marching procedures. 

The unique character of Matsuno’s scheme is that the crosswise (y) derivatives are 
formed independent of the sign of the crosswise velocity component. The crosswise deriva- 
tives are approximated by explicit, three-point central differencing at the previous step, 
which yields stability independent of the crossflow direction. Therefore, Matsuno’s finite- 
difference molecule does not depend on the crossflow direction. 

The zone of dependence principle requires 


u 

— >0 and 

u e 


hi Ax v 
hi Ay u 


(96) 


Matsuno’s finite-difference scheme is conditionally stable [16], and the stability condition 


gives the same constraint as that required by the zone of dependence principle. 


5.2.3 Accuracy 

The accuracy of the present procedure is established by comparing numerical results 
for several test problems with previously published results obtained by other investigators. 
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Case 1. Flat plate with attached cylinder 


Three-dimensional incompressible laminar flow past a flat plate with an attached cylin- 
der (Fig. 4) was calculated using Cartesian coordinates (hi = h^ = 1). This flow has 
been computed by several investigators. Fillo and Burbank [22] used the Crank-Nicolson 
method, Cebeci [23] used the Box scheme, and Iyer [8] used fourth-order scheme. For this 
flow, the inviscid velocity distribution is given by [23] 


„ f, , 2 y 2 -(x-x 0 ) 2 1 

• “l + 0 1 (x- *o) ! + !,’!’/ 

(97-a) 

.. _ OT/ -2 y(x-x o) 

2V ~» 0[(*_ 10 )> +y i]» 

(97.b) 


where ao is the cylinder radius, and xo is distance of the cylinder axis from the leading 
edge. To compare the numerical results with those obtained by other investigators, the 
conditions — 3050 cm/ sec, a$ = 6.1cm, and Xq = 45.7cm are chosen. The grid spacings 
are Ax — 0.61cm, Ay = 0.61cm, and Af = 0.2 with — 8.0. The results ((/ff)™) are 
shown in Table 1. The values presented in Ref. [23] have been multiplied by l/\/2 to 
properly account for differences in the transformation. The numerical results from the 
present method are in good agreement with those computed by Fillo and Burbank [22], 
by Cebeci [23], and by Iyer [8] as shown in Table 1. 


Table 1. Comparison of the values of ( f (( ) w 
y = 0 cm 


x(cm) 

Fillo and Burbank [22] 

Cebeci [23] 

Iyer [8] 

Present 


(C-N scheme) 

(Box Scheme) 


(Matsuno) 

0.00 

0.3321 

0.3319 

0.3321 

0.3323* 

2.44 

0.3292 

0.3289 

0.3289 

0.3293 

4.88 

0.3251 

0.3250 

0.3248 

0.3253 

7.32 

0.3199 

0.3198 

0.3195 

0.3202 
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9.76 

0.3130 

0.3130 

0.3126 

0.3133 

12.20 

0.3035 

0.3036 

0.3031 

0.3039 

14.64 

0.2903 

0.2907 

0.2900 

0.2908 

17.08 

0.2715 

0.2722 

0.2713 

0.2721 


y == 3.05 cm 


0.00 

0.3321 

0.3319 

0.3321 

0.3323 

2.44 

0.3292 

0.3289 

0.3290 

0.3294 

4.88 

0.3254 

0.3251 

0.3251 

0.3256 

7.32 

0.3203 

0.3202 

0.3200 

0.3206 

9.76 

0.3137 

0.3136 

0.3134 

0.3140 

12.20 

0.3048 

0.3047 

0.3045 

0.3052 

14.64 

0.2925 

0.2925 

0.2923 

0.2930 

17.08 

0.2751 

0.2751 

0.2752 

0.2758 


y — 6.10 cm 




0.00 

0.3321 

0.3319 

0.3321 

0.3323 

2.44 

0.3295 

0.3292 

0.3293 

0.3296 

4.88 

0.3260 

0.3257 

0.3257 

0.3262 

7.32 

0.3216 

0.3213 

0.3213 

0.3218 

9.76 

0.3159 

0.3156 

0.3156 

0.3162 

12.20 

0.3084 

0.3082 

0.3082 

0.3088 

14.64 

0.2985 

0.2985 

0.2984 

0.2990 

17.08 

0.2851 

0.2853 

0.2854 

0.2857 


where * is obtained using Blottner’s iterative method. 


28 



Case 2. Ellipsoid of Revolut ion 


The three-dimensional incompressible laminar flow over an ellipsoid of revolution with 
ellipticity ratio 4 (a = lm, b = 0.25m) was calculated for 0 and 6 degrees angle of attack 
using body-oriented coordinates and the analytical potential solution. This flow has been 
computed by Wang [24], Hirsh and Cebeci [25], and Cebeci and Su [26]. For this body, the 
metric coefficients can be obtained exactly, and the velocity components can be obtained 
analytically for the incompressible flow [10]: 


f l + (X/q-l)*(t 2 -l ) ) 1/2 
1 \ l-(X/q- l)* / 

h 2 = t\jl — JxJa — l ) 2 

u e — K»(ko(t) cos a cos /? — Vgo(t) sin a sin /? cos <f>) 
v e = Voo(Vgo(0 sinasin<£) 


(98. a) 
(98. b) 
(98.c) 
(98. d) 


where t = 6/a. Here (3 is the angle between the line tangent to the ellipse and the positive 
X axis; it is given by 


cos (3 — 

0< 0 


y/l -(*/»-!)» 

if X/a> 1, and (3 > 0 


if X/a < 1 


(99. a) 
(99. b) 


The parameters ko(t) and Vgo(0 are functions of t and are defined by 


Vo(0 = 


V 90 (t) = 


(1 - 1 2 ) 3 ' 2 

2V 0 {t) 

2V 0 {t) - 1 


(100. a) 
(100. b) 


The skin friction coefficients [C {x00 y/R^ 0 = where a = 1) as a function 

of X at an angle of attack zero degrees (axisymmetric flow) are shown in Fig. 5. The 
present numerical results were obtained using the following grid distributions: 40 steps of 
Ax=0.001 near the nose followed by Ax=0.02 downstream, and 41 grid points of Af =0.2. 
The present results are in very good agreement with the results of Hirsh and Cebeci. 
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However, the results of Wang are considerably different (higher) from the present results 
and also from the results of Hirsh and Cebeci. 

Table 2 shows the present results of the streamwise wall shear values, at an angle 

of attack 6 degrees at X = 0.5 (m). The present results were obtained using Ax=0.001 
for 40 steps near the nose, 0.02 downstream, Ay = 7r/36(kmax = 37), and Af=0.2. The 
present results and the Box scheme results computed by Cebeci and Su [26] are compared 
in Table 2. The results obtained by Cebeci and Su are multiplied by (u e /V 0 0 )~ 3 ^ 2 to cor- 
rectly account for the different definition of /' and transformation. The difference between 
using Box scheme and Characteristic Box is within 0.6 percent for Ax = 0.025. The dif- 
ference between the present result and the standard Box scheme is also within 0.6 percent. 


Table 2. Comparison of the values of /” at X = 0.5(m), a = 6° 


<f>(degree ) Standard Box [26] 

Characteristic Box [26] 

Present 

0 

0.6735 

0.6735 

0.6701 

20 

0.6665 

0.6676 

0.6626 

40 

0.6443 

0.6461 

0.6410 

60 

0.6090 

0.6115 

0.6064 

80 

0.5630 

0.5660 

0.5613 

100 

0.5103 

0.5134 

0.5096 

120 

0.4569 

0.4597 

0.4570 

140 

0.4116 

0.4138 

0.4121 

145 

0.4047 

0.4050 

0.4032 

150 

0.3952 

0.3969 

0.3956 

155 

0.3908 

0.3909 

0.3893 

160 

0.3837 

0.3852 

0.3843 

180 

0.3773 

0.3773 

0.3766 
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5.2.4 Efficiency 


Matsuno’s finite-difference method is fast and efficient compared with other current 
methods. The efficiency and speed of the Box scheme and Matsuno’s scheme can be 
directly compared. The Box scheme uses the block elimination method to solve the 6x6 
block tridiagonal system obtained from the momentum equations. The energy equation 
becomes a 2x2 block tridiagonal system when using the Box scheme. For Matsuno’s scheme, 
the momentum equations yield a 2x2 block tridiagonal system which can be efficiently 
solved using the Davis Modified Tridiagonal Algorithm, and the energy equation becomes 
a linear tridiagonal matrix form which can be solved by the Thomas algorithm. Another 
CPU advantage is that Matsuno’s finite-difference method is noniterative ( being only 
a predictor-corrector procedure) compared with the Box scheme, which is an iterative 
method because of the linearization. Although Matsuno’s finite-difference method requires 
smaller stepsizes (Az) near the stagnation point where the velocity gradients are high 
(which will be discussed in detail in Section 6.4), the computation time required to advance 
the solution to a given x station using Matsuno’s finite-difference method is substantially 
less than that required for the Box scheme. 

Matsuno’s finite-difference method is also highly vectorizable [16] compared with other 
schemes. The comparison of the CPU time between Matsuno’s and Box scheme can be 
found in Ref. [16]. According to this comparison, Matsuno’s scheme operating in the vec- 
tor mode requires approximately 1/50 of the CPU time required by the Box scheme to 
calculate the same number of grid points. For the cases presented in the present paper the 
CPU time per grid point was approximately 8xl0 -5 second on the CRAY-2 ( 8 seconds for 
110x31x31 grid ). 


5.3 Modification of Matsuno’s Finite Difference Method 

Matsuno’s finite-difference method uses explicit central differences for the crosswise 
derivative terms ( d jdy). Therefore if the solution at one of the side boundaries does not 
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exist then grid points are dropped as the solution is marched in x. In this case, two grid 
points are dropped (one at the predictor and the other at the corrector) for each step in the 
x direction (See Fig. 6(a)). Consequently, without modification, the Matsuno procedure 
can obtain the solution in only a small part of the whole flow field if the boundary-layer 
solution does not exist at any y grid point. 

The modified Matsuno’s finite-difference method uses second-order backward differ- 
ences for the crosswise derivative terms at the side boundary when the solution at one 
side (in the y-direction ) of the solution domain does not exist. The purpose of using 
this method is to continue the solution downstream while minimizing the number of lost 
solution stations in the x = const plane. This method can be used provided all the local 
streamlines are in the positive y-direction, i.e., the zone of dependence principle is satis- 
fied. When using this procedure, the boundary-layer solution can be obtained as far as 
the boundary-layer assumption is valid (see Fig. 6(b)). This procedure can be used even 
when the open type of separation occurs off the planes of symmetry. 

Figure 7 shows the difference molecule for the modified Matsuno’s finite-difference 
method. The modified finite-difference molecule is used only for the side boundary grid 
point. As shown in Fig. 6(b), the standard Matsuno’s finite-difference molecule is used for 
all interior points. 

The hypersonic flow over a cone with half cone angle of 10 degrees at 4 degrees angle 
of attack was chosen to validate the modified Matsuno method. The boundary-layer so- 
lution on the leeward line of symmetry ( <j> = tt) does not exist for these test conditions; 
consequently the modified Matsuno’s molecule will be used for y = y ma x- 
The flow conditions are the same as Tracy’s [27]: 

Moo = 7.95 

T too = 755.4° K 

P too = 1.7878x10 6 N/m 2 

T w /T t oo = 0.41 
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The inviscid solution was obtained using the Euler code developed by Manuel D. Salas 
(unpublished work) at the NASA Langley Research Center. The initial boundary-layer 
solution near the nose tip ( X = 0.02 m) from the windward line of symmetry to near 
the leeward line of symmetry was obtained using Blottner’s iterative method (Section 
5 . 1 ). Using these initial velocity and temperature profiles, the calculation was continued 
140 steps downstream to X = 2.39 m using the modified Matsuno procedure on the side 
boundary. This calculation was done using the body-oriented coordinate system. The step 
size (Ax) was small ( around 0 . 0007 ) near the nose tip, due to the zone of dependence 
requirements, and increased as the solution proceeded downstream to a maximum value 
0 . 2 . The number of grid points in the y and f directions are 31 (Ay = tt / 30 ), and 41 
(with Af = 0 . 2 ), respectively. The heat transfer ratio at the initial (J = 1, X = 0.02 m) 
and final ( I = 140 , X = 2.39 m) solution stations is presented in Fig. 8. The initial heat 
transfer at I = 1 ( X — 0.02 m ) is plotted as a line in Fig. 8; the heat transfer at I = 140 
( X = 2.39 m) is plotted as circles. The numerical results by Boericke [ 29 ], who obtained 
the similarity solution using Blottner’s iterative method based on the inviscid solution from 
Moretti [ 30 ], are also compared. The numerical results agree very well with each other as 
well as with the experimental data obtained by Tracy [27]. Consequently, the modified 
Matsuno procedure can be used to obtain accurate downstream solutions for those cases 
where the boundary-layer solution does not exist for a side boundary. 


33 



6. RESULTS AND DISCUSSION 


6.1 Solution Procedure 

The present method can be applied to any general fuselage-like configuration for com- 
pressible, perfect-gas flow. The calculations can be made for either body-oriented or 
streamline coordinates. A schematic flow chart of the procedure is presented in Fig. 9. 

Two geometry programs have been used for modeling the bodies studied in the present 
report: (1) the QUICK geometry program [31] and (2) the semi-analytic geometry program 
developed by Barger and Adams [32]. The QUICK geometry program was to used for the 
ellipsoid of revolution. Using the QUICK geometry program for a simple body, like an 
ellipsoid of revolution, is exact since it uses analytic functions for the arc and bodyline 
modeling. For the fuselage which was chosen for a test case, the semi-analytic geometry 
program developed by Barger and Adams [32] was used. For this nonanalytic body, the 
QUICK geometry program was not used due to the time it would have required to setup 
for the geometric model. The semi-analytic geometry program [32] was found to yield 
accurate and smooth modeling for the fuselage. 

The inviscid solutions used are (1) an analytic potential solution, (2) Euler code, and 
(3) the potential code developed by Hess [33]. Analytic potential solutions have been used 
for the incompressible flow over a flat plate with an attached cylinder and for the ellipsoid 
of revolution. In these cases, the analytically obtained inviscid solution and the metric 
coefficients can be given directly to the boundary-layer code. The analytical inviscid 
solution was used to investigate the accuracy of the finite-difference method in Section 
5.2.3. For the inviscid solution on the cone, Euler code developed by Manuel D. Salas 
(unpublished work) at the NASA Langley Research Center was used. The potential flow 
code developed by Hess [33] was used to obtain the inviscid flow field over the ellipsoid of 
revolution and the fuselage. 

A typical inviscid grid on the fuselage is shown in Figure 10. For the present study, 
most of the inviscid solutions from the Hess code were obtained using 54 grid points in 
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the X-direction and 37 grid points in the (^-direction. The Hess code gives the inviscid 
Cartesian velocity components and pressure coefficients at the centroids of each panel. A 
major problem with the panel method is the loss of accuracy in the nose region because of 
the singularity on the axis as X — * 0. Consequently, the boundary-layer calculations must 
be started slightly downstream from X = 0. 

Two computer programs were developed to calculate the boundary-layer coordinates 
(body-oriented and streamline ). These programs read the Cartesian inviscid velocity 
components and the pressure coefficients on the inviscid grid. Given x and y distribution, 
boundary-layer grid of the body surface is calculated using the method presented in Ap- 
pendix D.l (for the body-oriented coordinate system) or Appendix D.2 (for the streamline 
coordinate system). These programs calculate the following on the boundary-layer grid: 
u e , v e , s, cos 0, hi, hj, Cp. The velocity components, cos 0, and the pressure coefficients are 
interpolated from the inviscid grid onto the boundary-layer grid using bidirectional cubic 
splines with tension interpolation subroutine. 

For subsonic flows, the pressure on the body surface is not required as input because it 
can be calculated using the velocity components and the isentropic relationship with the 
freestream. When the pressure coefficients are not given on the boundary-layer grid for the 
subsonic flow, the three-dimensional boundary-layer code calculates the pressure using the 
isentropic relationship (see Appendix D.3 for detail). However, for the supersonic flows, 
the pressure coefficients on the boundary- layer grid must be specified because the pressure 
on the body surface is not related isentropicly to the undisturbed free stream. 

The boundary-layer calculation starts near the stagnation point or near the nose tip 
for the fuselage with the initial velocity and temperature profiles. The initial profiles 
are obtained at i = 1, and the boundary-layer calculation starts from i = 2 (See Fig. 
11). Each step (predictor or corrector step), the initial calculation starts at the windward 
line of symmetry (j — 1); the unknown points off the line of symmetry are solved for 
increasing values of j ( j — 2,3 ,..,jmax — 1); then the solution at the leeward line of 
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sy mm etry (j = jmax) is obtained. The present finite-difference method uses explicit 
central differences for the y-derivative terms. Consequently, the procedure does not sweep 
in the y-direction, as is required for most other procedures. Therefore, identical numerical 
results are obtained if the procedure is reversed, i.e., start at the leeward line of symmetry 
(j = jmax ) and solve for decreasing values of j (j = jmax — l,ymax — 2, ..,3, 2) and then 
the windward line of symmetry ( j — 1). 

The procedure described above is based on the assumption that the boundary-layer 
solutions on the whole surface (up to the leeward line of symmetry) exist. This assumption 
is generally valid before the separation line (closed or open type of separation) for a blunted 
nose fuselage. 

For a sharp nose fuselage, the initial solution near the nose tip (t = 1) is obtained 
using Blottner’s procedure [15] as far in the direction toward the leeward line of symmetry 
(j = 1, 2, ..) as the solution can be obtained. If the boundary-layer solution on the leeward 
line of symmetry exists, then the solution procedure for the sharp nose fuselage is the 
same as for the blunted nose fuselage. However, if the boundary-layer solution near the 
leeward line of symmetry does not exist, the modified Matsuno’s finite-difference method 
introduced in the Section.5.3 is used for the last point where the boundary-layer solution 
exists. 

In the present method, the coefficients, mi, m 2 , are determined numerically 

from the given velocity components (u e ,v e ), s, cos 9, and the metric coefficients (hi, hi). 
The coefficients, m x , m 2 , .., mu, are evaluated at the mid point (x,- +x / 2 ) except for the 
crosswise (y-) derivative terms, which are evaluated at the corrector step (x,+ x ). For ex- 
ample, (1 /h 1 )(du e /dx) is obtained using central differencing at the mid point (x, +x / 2 ), 
i.e., 

1 ^ due _ u < +1 ~ U 'e (101) 

h\ dx ds <Si+i iSj 

Similar derivatives, such as (l/hi)(dv e /dx), (l/h x )(d/i 2 cos 0/dx) are obtained in the same 
way as the above. However, the crosswise(y-) derivative terms, such as dujdy, dv e /dy, 
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dhi/dy, are obtained at the corrector step (x,+ i) assuming nonuniform spacing in the 


y-direction, i.e., 

du e _ (Ayj-i/Ayj^Uej+i — Uej) + (Ayy/Ay f --n)(u g j — u 
dy &yj + Ayy_ i 

The nonderivative properties, tt e , v t , s, h\, h 2 , p , and p are averaged 

t-th and i + 1-th step. 


hlzi L (102) 

from the values at 


6.2 Boundary-layer Parameters 


The skin friction coefficients are defined and calculated from: 

(pdu/dz) w _ 

/- " ±p e V> ~ Pe V> 

( pdv/dz) w _ 2/x UJ v 0O (aG/ac) u <(p/p e )«»(p e u t /M« 5 ) 1/2 

p*v? 


r - ^ 

C fy - — J 


\p.V* 


(103. a) 
(103. b) 


To compare the skin friction coefficients with the results obtained by other investigators, 
another definition of the skin friction coefficients is sometimes needed. The following 
definition is used for the incompressible flow over the ellipsoid of revolution: 

(pdu/ dz) w 2/X tt |U e (^-F , /3<r) tu (p/p e )m(PtUe//^c’S) ^ 


C/xoo — i 


r - 

/voo_ Sp^v? 


2PooV£ PooVl 

{pdv/dz) w _ 2/^V^p [dG/ds) w {p/pe)w(peUe/Hes) 1/2 

PooVl 


(104. a) 
(104. b) 


2 r oo r oo 

where (< BFjd f)«„ and (« dG/d $) w are evaluated by second order one-sided differences at the 
wall, i.e., 


(¥)- = 


__ (Aft + Aft) 2 F 2 -(Aft) 2 F 3 

dc ,w Aft(Aft + Aft) 2 - (Aft + Aft)(Aft) 2 

,<9G. _ (Aft + Aft) 2 G 2 — ( Aft) 2 G 3 

(^t)” " 


(105. a) 
(105. b) 


a? )w Aft(Aft + Aft) 2 - (Aft + Aft) (Aft) 2 

The skin friction coefficients presented in the present report are referenced to the body- 
oriented coordinate system. Results obtained in the streamline coordinate system have 
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been transformed to the body-oriented coordinate system as follows; 


{Cfz)b - (c% + Cj y ) ){ 2 {cos(fi + 7 ) - sin(/? + 7 ) cot 6 } (106.a) 

{Cfv)b = (C% + Cf v ))( 2 sin(/? + 7 ) csc (9 (I06.b) 


where /? = tan 1 {Cf v /Cf x ), t and 7 is the angle between the streamline coordinate line 
(y = const ) and the body-oriented coordinate line (y = const). 

Displacement thickness as presented in the present paper was not obtained from the 
displacement surface equation but, instead, from the following definition: 

6* = r< 1 _ fL) dz = rli- ( p ~ ) + (KooG,)2 + 2u ^ fg cos a] 1 /* 

Jo p e v/ Jo \ K pJ v e 

Heat transfer is calculated from: 



„ hrl^' F \ c p tLvj f P \ /Pe u e\ 1/2 

K ( dy )"~ p r ( / , e )“'( AieS ) ' ( 108 ) 

where (dT/d^) w is obtained from: 

_ [(Aft ) 2 ~ (Aft + Aft) 2 ]^ + (Aft + A ft) 2 T 2 - (Aft) 2 T 3 
d $ W Aft (Aft + Aft ) 2 - (Aft + Aft) (Aft ) 2 1 > 

Along the lines of symmetry, C fx , C fv , 6\ and q w are obtained by substituting G = 0 
in the equations above. 


6.3 Test Cases 


6.3.1 Hypersonic Cone With Mass Transfer 

The hypersonic flow over the sharp cone at 0 degrees angle of attack with mass transfer 
at the wall was selected as a test case with the following flow conditions (This flow condition 
is the same as that used in Ref. [ 8 ] and [28]): 

M 0 o = 7.4 

e c = 5° 


a = 0° 
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p too = 4.14xl0 6 Njm 2 
T too = 833 °K 
Tv = 316.65° 

From X = 0 to X = 0.096 m, there is no mass transfer at the wall. From X = 0.096 m, 
three types of mass transfer exist: (put). = 0; (p»). = -0.0901 N sec/ m? (wall section) ; 

and (pu>) w = —0.090117 N sec/ m s (wall injection). 

The inviscid solution was obtained using the Euler code developed by M. D. Salas, 
(unpublished) Figure 12 shows the skin friction coefficients on the cone with the mass 
transfer conditions listed above. These results were obtained using uniform grid spacing 
in the normal direction with Af = 0.2. From X = 0.096 to X = 0.1, very small stepsizes 
(Ax = 0.0002) were used to obtain a smooth skin friction coefficient when there is a mass 
transfer. The results are in good agreement with the other results [8] (not shown). Inci- 
dently, the solution obtained using a nonuniform (stretched) grid spacing in the normal 
direction (Af(l) - 0.01, A ${j + l)/A?(j) = 1.05, j=l,2,..,jmax-l) is also presented as 
small circles in Fig. 12. This result is in good agreement with the result using the uniform 
grid spacing in the f-direction (solid line). 

6.3.2 Supersonic Cone 

The supersonic flow over the sharp cone with a half cone angle of 5 degrees at an angle 
of attack 2.25 degrees was selected as a test case with the following flow conditions: 

Moo = 3.5 
0c = 5° 
a = 2.25° 

p too = 3.6x10 s lbf/ft 2 
T too = 540° 

T w = T a v 

The inviscid solution was obtained using the Euler code developed by M. D. Salas. The 


39 



initial boundary-layer solutions near the nose tip ( X = 0.02 ft) from the windward line 
of symmetry to near the leeward line of symmetry are obtained using Blottner’s iterative 
method (Section 5.1). The step size (Ax) is small ( around 0.0005) near the nose tip because 
of the zone of dependence requirements and was increased downstream. The number of grid 
points in the y and f directions are 31 (Ay = 7r/30) and 41 (with Af = 0.2), respectively. 

The skin friction coefficients at X = 0.124,0.208,0.348, and 0.582 ft are plotted in 
Fig. 13. At this angle of attack, the boundary-layer solution along the leeward line of 
symmetry does not exist. However, using the modified procedure presented in section 
5.3, the solution could be obtained downstream. The boundary-layer thickness and the 
displacement thickness at the same locations are plotted in Figs. 14 and 15, respectively. 
The streamwise and crosswise velocity and the temperature profiles at X = 0.582 ft are 
plotted in Figs. 16 through 18. 

6.3.3 Ellipsoid of Revolution 

An ellipsoid of revolution having a four to one ratio of major to minor axis (a = 1 m, 
b = l/4m) was selected as a test case, and the boundary-layer solutions were obtained 
for incompressible flow at a = 6°. This particular case was selected because (1) its geom- 
etry is analytic, (2) an exact potential solution exists for the inviscid flow field, and (3) 
numerical results have been previously published (Reference [34]). Two approaches were 
taken for each of the coordinate systems: (1) analytic grid generation, analytic metrics, 
and the analytical inviscid potential solution ;(2) numerical grid generation, numerically 
calculated metrics, and the inviscid flow field from the panel method of Reference [33]. The 
axisymmetric analogue [35] was also used to obtain approximate results for comparison. 

Skin friction coefficients (here, C fxoo y/ReZ = where a = 1) and the 

displacement thickness results for = 1 m/sec are presented in Figures 19 and 20, re- 
spectively. The skin friction results presented in Figure 19 are in the body-oriented coor- 
dinate system; i.e., results obtained in the streamline coordinate system were transformed 
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to body-oriented coordinate system using Eq. (105) . When the exact metrics and inviscid 
potential solutions were used, the difference between the results obtained using the body- 
oriented and streamline coordinate systems was less than 0.25 percent. The agreement 
between the numerical results obtained using the exact inviscid potential solution and 
analytic geometry and the more general numerical approach, i.e., numerically generated 
inviscid solution and coordinate system, is excellent for both boundary-layer coordinate 
systems. In addition, comparisons are made with the results presented in Reference [34]. 
The axisymmetric analogue results have the correct trend as compared with the three- 
dimensional results, except near the three-dimensional separation line, and are generally 
within ±5 percent of the three-dimensional values away from the separation line. At this 
angle of attack, the three-dimensional separation line begins approximately at X = 1.32 m, 
<f> = 110° (See Cebeci and Su [26]). 

6.3.4 General Aviation Fuselage 

A low speed, general aviation aircraft fuselage was selected as a test case with nonana- 
lytic geometry. The particular configuration has served as a flight test vehicle for transition 
prediction procedures (see Refs. [36] and [37].). This case is particularly interesting in 
that the crossflow is into the plane of symmetry: < 0 as <£ — ►0,t^>0as<£— * tt. 

Consequently, standard marching procedures can not be used to advance from the solution 
obtained on either line of the symmetry plane (<f> = 0, <f> = 7r) into the three-dimensional 
region (0 < <f> < 7 r), because any attempt to do so would violate the zone of dependence 
principle. 

Numerical results were obtained for a Mach number and unit Reynolds number of 0.3 
and 7xl0 6 m“\ respectively for 0° and 3° angles of attack for an adiabatic wall. A photo- 
graph of the aircraft is presented in Fig. 4 of Reference [36]. A typical panel distribution 
used to obtain the inviscid solution is presented in Figure 10. The boundary-layer grid for 
the body-oriented coordinate system is shown in Figure 21, and two boundary-layer grids 
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for the streamline coordinates (31 streamlines) are shown in Figure 22. 

(1) Zero Degrees Angle of Attack 

A side view of the fuselage forebody is presented in Figure 23 showing the maximum 
pressure line and the crossflow velocity regions. As previously discussed, because the lines 
of symmetry = 0, and <f> = 7 r) are inflow lines, standard marching procedures can 
not be used to advance the solution from the lines of symmetry into the interior region 
(x = const; j — 2,3 ,..,jmax — 1). However, Matsuno’s method is independent of the sign 
of the crossflow velocity and can be used. 

Numerical results are presented in Figure 24 through 26 for a — 0°. The agreement 
between the skin friction coefficients obtained in the two coordinate systems is excellent 
over the entire surface and the difference is within one percent. The axisymmetric analogue 
results have the same general trend as the three-dimensional boundary-layer results, except 
that they yield larger values of C/ x near the plane of symmetry and fail to predict the Of* 
trend along the side of the fuselage (<f> « tt/2). 

Boundary-layer thickness and displacement thickness results exhibit similar trends in 
agreement between the results obtained in the two coordinate system and the axisymmet- 
ric analogue. 

(2) Three Degrees Angle of Attack 

At this angle of attack, the flow field has two relative maxima pressure lines for 0 < 
<f> < 7r with multiple changes in the sign of the crossflow velocity (see Figure 27.). Figure 
28 shows the values of v e /Voo as a function of <j> from the inviscid solution at X = 0.6, 0.9, 
1.2, and 1.5. The sign of v e changes three times as <f> increases from 0 to n at X = 1.5. 
Numerical results for a = 3° are presented in Figures 29 through 35. Example streamwise 
and crossflow velocity profiles are presented in Figs. 29 and 30. Temperature profiles are 
presented in Figure 31. 

Skin friction coefficients, boundary- layer thickness, and the displacement thickness re- 
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suits are presented in Figures 32 through 34, respectively. The agreement between the 
numerical results in the body-oriented and the streamline coordinate system is excellent 
over the entire surface. The axisymmetric analogue results are not acceptable at this an- 
gle of attack in that they are in error in both magnitude and trend. For example, the 
axisymmetric analogue results for Cj x are on the order of 10 percent greater than the 
three-dimensional values as the windward plane of symmetry is approached and fail to 
predict the C j x trend for <f> « 7r/2. 

Three-dimensional results for the velocity profiles in the streamline direction are com- 
pared with the axisymmetric analogue results in Figure 35 for X = 1.3m. The axisymmet- 
ric analogue results yield larger values of d(uju e ) w /dz than the three-dimensional results 
as well as larger values of uju t across the boundary layer at these points. 

During the fuselage study questions were raised concerning the behavior of the numer- 
ical results in the neighborhood of the plane of symmetry. It can be seen that the results 
from the streamline coordinate system were not smooth; i.e., slight oscillations occured 
near the plane of symmetry. The streamline coordinate system using 31 streamlines for 
the fuselage is presented in Figure 22. The streamlines are initially uniformly distributed 
in the crossflow plane at the initial station (x « 0); however, they tend to converge to- 
ward the symmetry lines (<j) = 0; (f> = 7r) as x increases. Consequently, as x increases, 
the streamline distribution becomes highly nonuniform with a dense packing of the grid 
lines in the neighborhood of the symmetry lines and a sparse distribution in the region 
tt/3 < (j> < 5tt/ 6. Using this streamline distribution it was not possible to obtain the 
correct skin friction coefficient behavior numerically near <f> = 7 t/ 2. This problem was 
not present in the body-oriented coordinate system, because the boundary-layer grid re- 
mained nearly uniform for each x location. In order to obtain the correct behavior of the 
boundary-layer in the (f> « 7r/2 region it was necessary to use 91 streamlines. Although 
the use of 91 streamlines corrected the problem in the region of <f> & vr/2, it created a new 
problem, oscillatory values of the boundary- layer parameters, in the neighborhood of the 
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symmetry lines. This was found to be caused by the combined effect of grid density and 
the oscillation of the metric coefficient, / 12 , in the y-direction. Using 31 streamlines, this 
oscillation was not present. 

The term dKi/dy assumes a dominant role along the symmetry lines. The term is 
directly proportional to d 2 u e /dy 2 because the metric coefficient, hi, is defined by Eq. (23) 
for the streamline coordinate system. The inviscid results are first obtained at the cen- 
troids of the panels. Using this relatively coarse distribution of centroids, it is difficult 
to obtain accurate values of d 2 u t jdy l on the densely packed boundary-layer grid because 
due/dy is zero at <f> = 0 and <f> = n. Special attention must be paid to this term when 
using the streamline coordinates and when the grid points are concentrated near the lines 
of symmetry. 


6.4 Restricton on Grids 

The present method is developed in such a way that one can use nonuniform grid 
spacing in the streamwise (x), crosswise (y), and in the normal (f) direction. The grid 
spacing in the x-direction near the stagnation point where the velocity gradients are large, 
has to be small not to yield oscillatory boundary-layer parameters. Figure 36 shows the 
computed skin friction coefficient variation with X for the incompressible flow over the 
ellipsoid of revolution (a = lm,6 = 1/4 m) at zero degrees angle of attack. The analytic 
inviscid solution as given in Section 5.2.3 is used to obtain the results. Step size (Ax) values 
of 0.001, 0.002, and 0.005 were used for the first 40 steps; thereafter, Ax was set to 0.02. 
As can be seen in this figure, the computed skin friction coefficient near the stagnation 
point oscillates increasingly about the correct solution as the step size is increased. In this 
figure, the correct solution can be assumed to be the solution of using Ax of 0.001 near 
the stagnation point (solid line). It is also apparent that once the solution approaches the 
correct solution the oscilation vanishes even with abrupt increase in Ax to 0.02. 

The zone of dependence principle requirement (Eq.(96)) also restricts Ax. Because 
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the coordinates are calculated prior to the boundary-layer calculation (see Fig. 9), the 
stepsizes are determined before the stability condition is checked. The stability is checked 
as the boundary-layer solution is obtained, and the calculation is continued as far as the 
stability condition is satisfied and is stopped when it is violated. Therefore, when the zone 
of dependence principle is violated with the present Ax distribution, the stepsizes after 
that point must be obtained by a trial-and-error method. However, the step where the 
stability condition is violated must be very close to the separation line because the wall 
limiting streamline changes its direction very rapidly near the three-dimensional separation 
line. It is to be noted that the stepsize near the nose is severely restricted because of the 
zone of dependence principle. 

The grid distribution in the f direction used to obtain most of the results presented was 
uniform (for most cases, Af=0.2). However, a nonuniform grid spacing in the f direction 
can be used as well. The boundary-layer grid in the y direction can be given arbitrarily. 
However, the grid points in the y direction used in this study are uniformly distributed 
(Ay = n/(kmax— 1)) regardless of the coordinate system used (body-oriented or streamline 
coordinates). The actual distance between two grid points with the same x is h 2 Ay. In the 
streamline coordinates, y remains the same along each streamline; therefore, Ay between 
two streamlines remains the same even downstream. The metric coefficient, h 2 , is a direct 
function of the streamline divergence and varies over a large range as x increases in the 
streamline coordinate system. Nonuniform spacing downstream on the general fuselage 
when using the streamline coordinates is due to the variation of metric coefficient h 2 , not 
because of the nonuniform distribution of Ay. 
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7. CONCLUDING REMARKS 


The three-dimensional, compressible, laminar boundary-layer equations have been nu- 
merically solved for several configurations at angle of attack. The finite-difference pro- 
cedure used to solve the governing equations is second order accurate. The procedure is 
independent of the sign of the crossflow velocity component and is the best method known 
for configurations having crossflow reversal regions. The crossflow velocity direction for a 
general aviation fuselage was into the plane of symmetry on both the most windward and 
leeward surfaces for the angles of attack considered (0° and 3°). Consequently, standard 
solution procedures, which march around the body using the plane of symmetry as an 
initial data plane, could not be used to solve this fuselage test case. However, no numer- 
ical problems were encountered using the finite-difference procedure used in the present 
analysis. This was true even at three degrees angle of attack, where the crossflow veloc- 
ity component reversed direction as many as three times in the region bounded by the 
windward and leeward symmetry planes. 

Numerical solutions for the fuselage-type configurations were obtained using two boundary- 
layer coordinate systems: (1) a body-oriented coordinate system and (2) a streamline co- 
ordinate system. The agreement between the boundary-layer parameters obtained in the 
two coordinate systems was excellent over the entire fuselage surface. The boundary-layer 
grid (body-oriented or streamline) can be totally independent of the inviscid grid, i.e., the 
number of grid points in the x and y direction of the boundary-layer grid do not have to 
be the same as those of the grid used to obtain the inviscid solution. 

Based on the experience of using the two different coordinate systems on a general 
fuselage, the following conclusions can be made: (l)the generation of the streamline coor- 
dinates requires more effort than the body-oriented coordinates; (2) it is difficult, if not 
impossible, to control the boundary-layer grid distribution using the streamline coordi- 
nates; (3) the numerical results obtained using the body-oriented coordinates with only 31 
grid points in the y-direction are better than those using the streamline coordinates with 
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91 streamlines. Therefore, the body-oriented coordinate system is generally better than 
the streamline coordinate system for boundary-layer calculations on a general fuselage, 
provided the geometry singularity at the nose point in the body-oriented cooedinates can 
be avoided. Note, however, that the streamline coordinate results are excellent and may be 
preferable when output along streamlines is needed, as in transition prediction procedures. 

Excellent agreement of the boundary-layer solutions using two different coordinate sys- 
tems strongly validates this bound ary- layer method and the application software. This 
boundary-layer method is robust, fast, and can be applied to any type of fuselage (either 
with a blunted nose or sharp nose) which has a symmetry plane. However, further devel- 
opment is needed to add additional capabilities, such as viscous-inviscid interaction, real 
gas effects for hypersonic flows, and turbulence closure. A user’s manual with a detailed 
description of the computer programs and input is presented in Volume II. 
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Appendix A. Block Tridiagonal Matrix Algorithm 


The two vector equations to be solved are 

h k = h k . x + ^ ~{H k + H k _ 0 (Al.a) 

—A k H k -i + B k H k — C k H k+ 1 + a k h k = D k (Al.b) 

where A k , B k ,C k , and a k are 2x2 matrices, II k , h x , and D k are vectors. These equations 
are solved using the Davis Modified Tridiagonal Algorithm (See Appendix A). Introduce 
E k , e k and d k such that 

H k = E k H k - 1 + e k h k -i + d k (A2) 

where E k and e k are 2x2 matrices, d k is a vector. 

Using Eq. (A2), Eq. (Al.b) becomes 

— AjfciZjfc_i + B k H k — C k E k+ iH k — C k e k+ ih k — C k d k+ 1 + a k h k = D k (A3) 

Define 

R k = a k — C k e k+ 1 (A4) 

Then, Eq. (A3) may be written as 

— A k H k -i + ( B k — C k E k+ i)H k + R k h k = D k + C k d k+ i (A5) 

Substituting Eq. (Al.a) into Eq. (A5) gives 
(~A k + + ( B k - C k E k+ i + -R k )H k + R k h k _ x — D k — C k d k+1 = 0 (A6) 

Next, define 

p k — B k — C k E k+ 1 4 -R k (A7) 

Solving Eq. (A6) for H kf 

Hk = P* X (A* - R k )H k - 1 — p k 1 R k h k -i +p k l [D k + C k d k+X ) (A8) 
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Equating Eqs. (A2) and (A8) term by term yields 


(A9.a) 

(A9.b) 

(A9.c) 


e* = -P k lR k 

rr> -1a i ^Ck—1 

Ek = Pk fed — — ek 

dk = Pk'iDk + C k d k+ i) 


The boundary condition at the edge of the boundary-layer (k=kmax) is 


H 


kmax 


1.0 

v e /V ref 


This provides the conditions 


d 


kmax 


1.0 

v./V„ f 


ekmax = Ekmzz = 0 (A10) 

The parameters of Eq. (A9) are first determined for decreasing values of k (kmax-1, kmax- 
2,. ..,2) beginning at the edge of the boundary-layer. Then Eqs. (Al.b) and (A2) are solved 
for increasing values of k (k=2,3,..., kmax) using the boundary conditions at the wall, 

Hi = hi = 0 (All) 
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Appendix B. Initial Profiles Near Stagnation Piont 


Here, to avoid confusion, the rectangular coordinates, which have their origin at the 
stagnation point are described by x*, y*, and z* (See Figs. 37 and 38.) The corresponding 
velocity components are u*, v* , w*, and the inviscid velocity components are u*, v*, w*, 
respectively. The relations of this rectangular coordinate system ( x *, y* , z*) with the 
rectangular coordinate system (x', y', z ') which has its origin at the body nose ( X = 0), 
are 


x* = cos0 r (2' — z ' s ) — sin0 r (x' — x \ ) (Bl.a) 

y* = y' (Bi.b) 

z* = —cos0 r (x' — x' t ) — sin 0 r (z' — z' s ) (Bl.c) 

where 0 r is the angle between the two coordinate systems, and x' t and z' t are the x' and z' 
coordinates of the stagnation point, respectively. 

The stagnation point solutions are denoted with the subscript s, i.e., 

fl = u*/u* e (B2.a) 

9 = v*/v* e (B2.b) 

and are obtained using the Blottner’s iterative method. 

B.l Body-Oriented Coordinates 

The body-oriented coordinates (x, t/) which have their origin at X = 0 and the rect- 
angular coordinates ( x*,y *) which have their origin at the stagnation point are shown in 
Fig. 37(b). It is assumed that the stagnation point is very close to X = 0, which is true 
for a small angle of attack regardless of the shape of the nose. 

The inviscid velocity components u t and v e at a point P are 


u e = -u* cos y + v* e sin y 

v e = u* sin y + v* cos y 
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(B3.a) 

(B3.b) 



The boundary-layer velocity components at the point P are 


(B4.a) 

(B4.b) 


U = - u* /' cos y + v* e g' a sin y 
v = u*/' sin y + v* e g\ cos y 

The initial velocity profiles (off the lines of symmetry) in the body-oriented coordinate 

system are obtained from the following equations. 

, _ — u*f a cos y + v* e g', siny ^ cos y + By g a siny 

f =u/u e - _ u * cosy + v * s i ny —Ax* cos y + By* sin y 

, u! f' t sin y + v*g' t cos y A c*/i sin y + ffyV, cos y fB 5 b) 

9 = v/V„ jT V., 

where A and B are velocity gradients at the stagnation point (defined in Eq. (69)) in the 
x* and y* directions, respectively. 

On the lines of symmetry, 

r = u/u, = r, ( B6 ) 

Even though v is zero along the lines of symmetry, dv/dy is not zero and can be obtained 
as follows (y* = siny = 0 ): 

g' = VyjVoo = ^-(Ax* f' siny + By* g' cosy) /V^ (B7) 

oy 

or 

g ' = vjVoo = (Ax* f' cosy + B{dy*/dy)g'cosy)/V 00 (B 8 ) 


Note that 


cos y = 1, dy* jdy - h 2 along the windward line of symmetry (B9.a) 
and cos y = — 1, dy*/dy = -h 2 along the leeward line of symmetry (B9.b) 


Finally, 

g' = vJV m = (~A\x'\r, + Bh 2 g',)/V„ (BIO) 
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B.2 Streamline Coordinates 


The streamline coordinates ( x,y ) and the rectangular coordinates (x*,y*) which both 
have their origins at the stagnation point are shown in Fig. 38. 

The inviscid velocity u e at a point P is 


u e = u* sin /? + v* cos / ? 


(BH) 


The boundary-layer velocity components at point P are 

u = u lf, sin/? + v* e g' t cos /? 
v = Kf's cos /? ~ v* e g' a sin/? 


(B12.a) 

(B12.b) 


where /? is the angle between the streamline direction x and the y* direction and is given 
by 

* A * * 

X 


. a K Ax * 
tan 8 — — = = 

v* t By * C*y* 

Substituting (3 into Eqs. (Bll) and (B12) gives 

r = u/«, = f l t 

' 1 + {By-/Ax’Y 

a’ = „IV - 

' " v„([Ax-y + (Bz-yyn 

On the lines of symmetry, 

/' = u/u e - /' 

Along the lines of symmetry, dv/dy can be obtained from 

dv _ dv dx* dv dy* 

dy dx* dy dy* dy 

Near the lines of symmetry the following relations are valid: 

dx* 

- F - = h 2 cos/? 

dy 

dy* 

= —hi sin/? 
dy 


(B 13) 


(Bl4.a) 

(B14.b) 

(B15) 


(B16) 
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(B17.a) 

(Bl7.b) 


0 = 0; consequently, the first term of the right hand 


Along the lines of symmetry, cos 
side of Eq. (B16) vanishes to yield 

— = — h 2 sm (3-^-^{u* e f' s cos (3 — v*<j£sin/3) (B18) 

dy °y 

Expanding Eq. (B18) and substituting cos /? = 0, sin 2 /? = 1, dv*Jdy* = B,v e — By - 0, 
dpjdy* = -C*/x\ and u\ = Ax* gives 


g' = Vy/Voo = h 2 B(g' t — f',)/V oo 


(B19) 


The total enthalpy profile at point P is given as the stagnation point total enthalphy 
profile. 

An additional factor that must be considered is the difference between f in the stagna- 
tion point and main coordinate systems (body-oriented or streamline coordinate systems). 
This is because of the difference in definition of u e and s between the stagnation point 
coordinates and the main coordinates. To have a desired f distribution for the main co- 
ordinate system, interpolation could be used in obtaining the values at a corresponding 
f from the stagnation point solution. However, this procedure is not used in the present 
computer program because the difference of ujs is negligible when the angle of attack is 
small, as considered in the present report. 

The velocity gradients at the stagnation point (A, B ) and the location of the stagnation 
point (x',z;) and 0 r are calculated in the streamline coordinate program. To obtain good 
approximate initial profiles near the stagnation point using this procedure, the velocity 
gradients, location of the stagnation point, and 6 r must be obtained accurately. However, 
the velocity gradients are difficult to calculate accurately, especially when using a numerical 
inviscid solution, because the inviscid solution near the stagnation point changes rapidly. 
Accuracy was tested for the ellipsoid of revolution on which these values can be obtained 
exactly. These values could not be obtained accurately from the Hess panel method inviscid 
solution. The inaccuracy was caused mainly by the singularity of this inviscid method near 
X = 0. Thus the above procedure for calculating the initial profiles is used in the computer 
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program for the case when these values can be obtained accurately. 

However, when the angle of attack is small, as considered in the present study, the 
ratio of the two velocity gradients (B/A = C*) is close to unity and the stagnation point 
is close to X = 0. We can use the axisymmetric stagnation point solution obtained by 
using C* = 1 when any of the velocity gradients, location of the stagnation point, or 0 r is 
not easy to obtain or is not sufficiently accurate. In this case, the initial profiles near the 
stagnation point can be obtained easily. The difference of the boundary-layer solutions 
between cases using different values of C* vanishes within 5 downstream steps. 
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Appendix G. Initial Profiles Near Nose Tip for Sharp Nose Body 

The initial velocity and temperature profiles near the nose tip of the general sharp nose 
body are obtained based on the body-oriented coordinate system. Consequently, if the 
boundary-layer solution is sought using the body-oriented coordinates for the whole flow 
field on the sharp nose body, the initial profiles are used as they are obtained at t = 1. 

However, when the boundary-layer solution is to be obtained using the streamline coor- 
dinate system, initial velocity profiles based on the streamline coordinate system must be 
obtained, see Fig. 39. Assume the velocity profiles based on the body-oriented coordinates 
near the nose tip (obtained using the Blottner’s iterative method in the present report) 
are given as: 

fl = («/«*)» ( C1 - b ) 

g[ = (u/Voo)* ( Clc ) 


where subscript b denotes the body-oriented coordinates. 
Then, from Fig. 39, 


fat ~~ { U / U c)st ~ 


9',, = wvy.. = 


(,/u 2 + t> 2 )tcos/j 

{\fui + « e 2 )t 

-(\/u 2 + v 2 )t sin/? 


(C2.d) 


(C2.e) 


where the subscript st denotes the streamline coordinates. Here, (3 is the angle between 

the inviscid streamline and the local streamline and is given by 

(uu e + v V e )b 


cos /? = 


Substituting (3 into Eq. (C2) gives 


(\/u 2 + v 2 )i,(Vu 2 + u 2 ) t 


/J t = 


9 st 


(ufkfl + Vco(^) t g; 

(u 2 + U 2 )ft 

Voo{u e ) b g' b - 


(C3) 

(C4.g) 

(C4.h) 


Vooiy/ul + vl), 

In deriving Eqs. (C2) through (C4), the body-oriented coordinates ( x and y) near the nose 
tip are assumed to be orthogonal. 
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Appendix D. Boundary-Layer Edge Conditions 
D.l Body-Oriented Coordinates 

The body-oriented coordinate system employed in the present report is constructed by 
pure cross-sectional cuts. One coordinate is the line of intersection of the body surface and 
a X = const plane, the other coordinate on the surface is the line of intersection of the 
body-surface and a meridional plane (<^=const plane). This coordinate system ( x and y) 
is in general nonorthogonal for fuselage shapes. In these coordinates, x is measured along 
the axis and has the same value as X. Also, y has the same value as <j>. However, the 
directions of x and y are different from X and <f>, respectively. See Fig. 1 for the definitions 
of this coordinate system. Figure 21 shows the body-oriented boundary-layer grid on the 
fuselage shape studied in the present report. 

The inviscid velocity components in the rectangular coordinates are given at the cen- 
troid of panel P (X p ,<^ p ,r p ) from the inviscid code; see Fig. 40 . Using the geometry 
program, the point Q (X q ,4> q ,r q ), which lies within a small distance (6X) from the point 
P along x-direction, and the point R (X r ,<j> r ,r r ), which lies with a small angle {8<j>) from 
the point P along the positive y-direction, can be obtained. Here, 6X and 6<f> can be 
chosen arbitrarily small (typically 0.01). 

Now, from the definition of the rectangular (x', y', z') and cylindrical (X, r, <f>) coordi- 
nate systems, 


II 

y'p 

= r p sin (j) p 

z' P = -r p sin 4> p 

(Dl.a) 


y' q 

= r q sin 4> q 

z' q = -r, sin 4> q 

(Dl.b) 

x' r = X r 

y r 

= r r sin 4 *>r 

z* T = —r T sin 4> r 

(Dl.c) 


where X q X p + 6X , <f> q — <f> pt X r — X p , and <f> T = <f> p -\- 6(f>. 
Then, cos 6 is calculated from the following equation: 


K ~ - *;) + « - yp)(yr - £) + (4 - 0(4 - <) 

VK - z'r) 2 + (y; - 4) 2 + R - M - yr) 2 + R - kT 2 


(D2) 
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The angle between the streamline with the x-direction can be obtained by: 


cos 7 


= (u,./^)K - 4 ) + - yi.) + (D3j 

(V./^VK - x’ p Y + (sj - y- p ) ! + (*; - z’ p Y 


where V e is the inviscid total velocity, i.e. , 

K = yf‘ u[> + u®, + uf, ( D4 ) 

The inviscid velocity components in the body-oriented coordinates on the centriods are 
obtained from the equations below. 

Ue/Voo = (Vi/V'oo) sinqfcsc# (D5.a) 

v e /^oo = -(V e /V 00 )sin7 cotd+ (V' e /V 00 )cos7 (D5.b) 


After the inviscid velocity components and cosd are obtained on the centroids of panels, 
u e is extrapolated along the lines of symmetry. On these lines, cos 6 and v e are equated 
to zero. The velocity components, pressure coefficients (which is not necessary for the 
subsonic flow), and cos 0 on the boundary-layer grid are interpolated using the bidirectional 
cubic spline with tension program. The first derivatives of the velocity components, such as 
dug/dx , dug/dy, dv e /dx, are smooth and continuous using this interpolation subroutine. 
This subroutine must not be substituted by the Lagrangian interpolation subroutine. 

The metric coefficients are calculated as below: 


hi = 
h-i = 


ds ds 
dx dX 


IS* + (*£)* + (**)* 

{ dy ] K dy ] + W 


(D6.a) 

(D6.b) 


where dx'/dy , dy'/dy, and dz'/dy are obtained by central differences. 

For the case of the incompressible flow over the ellipsoid of revolution, the exact metrics 


and the analytical inviscid solution can be obtained in a closed form; see Section 5.2.3. 
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D.2 Streamline Coordinate System 


To calculate the inviscid streamlines, the method developed by Hamilton et al. [38] is 
used. In the present report, a spherical coordinate system ( R , 0, <fi) with the origin at 
X = X otp is employed (Fig. 41) for calculating the inviscid streamlines. The axial location 
X osp can have arbitrary value, but X otp = 1.0 is used. 

The velocity components in this spherical coordinate system can be obtained from the 
velocity components in the rectangular coordinates ( x ', y', z') by the following relationships. 

ur/Vco = {(*' - X 0ip ){u t ,/V 0 o) + y'(v/Voo) + 2 'K'/Ko)} /R (D7.a) 

«©/^oo = {y'{ x> — X oap ){u y > /V^) + z'(x' — X osp )(u z > /Voo) — r 2 (u I »/V 00 )| jrR (D7.b) 
u*/Voo = {z^Uyi/V oo) - 2 /(u*(/Foo)} /r 

where 

r = \f7 2 + 2 ' 2 

R = ^(x' - X 0ip ) 2 + + ** 

Now, let 

— = A 

be the derivative along an inviscid streamline on the surface. 

Then, from Eq. (50) and Ref. [38], 

— = h =* 22 . 

Dx 1 u e 

DX _ cos ©(ufl/Voo) - sinO^e/Voo) 

Dx K/Koo) 2 

-D<^> _ u^jVoo 

Dx r(u e /V oo) 2 

To establish the initial location of streamlines for the blunt nosed body, locate 0 4 (0 at 
the stagnation point) as the point where <f> = 0 and u e = 0. Then, draw a cone which has 
an angle e with respect to the line connecting the stagnation point and the origin of the 


(D7.c) 

(D8.a) 

(D8.b) 

(D9) 

(DIO. a) 
(DlO.b) 
(DlO.c) 
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spherical coordinates; see Fig. 41 . The initial location of streamlines on the intersection 
of the e cone and the body surface can be obtained from the following equations: 

= 7 r /2 + sin -1 | (e cos y cos 0, — Vl - e 2 sin0 3 )/Gl} (Dll. a) 

where 

^ /2 

G1 = | (e cosy cos 0, - Vi - e 2 sin©,) 2 4- sin 2 e sin 2 y} (Dll.b) 

cos 0 e = sin e cos y sin 0, + cos e cos 0, (D 1 1 .c) 

Equations (Dll) locate the coordinates ( X t and <f > t ) for the initial position of streamlines 
near the stagnation point. 

It should be noted that the initial locations of streamlines obtained as described above 
are not on the same x location, i.e., the initial x-direction and y-direction are not orthogonal 
to each other. To generate the orthogonal streamline coordinates, the initial location of 
streamlines has to be readjusted using integrations along the streamlines. 

For the sharp nose fuselage, the initial locations of streamlines are at the same small 
X. As for the blunted nose body, the initial locations of the streamlines ( i.e., on the same 
X) are not on the same x. To generate the orthogonal streamline coordinates, the initial 

locations of the streamlines have to be readjusted. 

The inviscid total velocity (u e ) and three velocity components (u R ,u^,u 0 ) are calcu- 
lated on the centroids of panels using Eq. (D7). The inviscid velocity components along 
the lines of symmetry are then extrapolated using the appropriate condition along these 
lines (u e ,u R ,u 0 ; symmetry condition). The inviscid velocity components and the pres- 
sure coefficients on the whole surface are obtained using the same interpolation program 
(bidirectional cubic spline with tension ) as used for the body-oriented coordinates. 

The fourth-order Runge-Kutta method is used to obtain the streamlines, i.e., to inte- 
grate Eq. (DIO). The metric coefficient, hi, is defined as Eq. (23) in this coordinates. The 
metric coefficient, h 2 , is obtained in the same way as for the body-oriented coordinates, 
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i.e., 


h<i — 




(D12) 


D.3 Temperature and Pressure 

For both body-oriented and streamline coordinate systems, the temperature at the edge 
of the boundary-layer can be obtained using the total inviscid velocity, i.e., 

T./T m = 1 + - (^)’| (D13) 

Equation (D13) is derived from the inviscid energy equation and is valid for all speed 
regimes. Temperature at the boundary-layer edge is calculated using the above equation 
in the boundary-layer code. Therefore, temperature at the edge of the boundary-layer is 
not a required input. 

For subsonic, shock-free flow, the pressure can be obtained from the isentropic relation 
with the free stream, i.e., 

p./p~ = (D14) 

oo 

Pressure at the edge of the boundary-layer is calculated using the above equation with T t 
obtained from Eq. (D13) in the boundary-layer code. Consequently, pressure at the edge 
of the boundary-layer is not a required input for subsonic flow. However, for supersonic 
flow when there is a shock wave present between the body and the free stream the pressure 
must be given as an input to the boundary-layer code, because the pressure on the body 
surface can not be obtained using the isentropic relationship with the undisturbed free 
stream. 
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Figure 1. Body-Oriented Coordinate System 
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Figure 3. Matsuno’s Finite-Difference Molecule 
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Figure 4. Flat Plan; with Attached Cylinder 
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Fig. 5 Skin Friction Coefficients (Ellipsoid of Revolution, = 1 m/sec, a = 0°) 
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Figure 6. Sketch for Modified Malsuno’s Procedure 
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Figure 7. Modified Matsuno’s Finite-Difference Molecule 






Figure 9. Flow Chart 
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Figure 10. Typical Inviscid Grids on a Fuselage 


75 




(b) Sharp Nose Fuselage 


Figure 11. Marching Procedure for Boundary-Layer Solution 
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Figure 13. Skin Fric'.ion Coefficients (Cone, — 3.5, a — 2.25°) 


78 



6 (ft) 


0.005 


0.004 


O 

% 


0.003 


0.002 


0.001 


0.000 


X (ft) 

0.124 

0.208 

0.348 


□ 0.582 


i o o o d o o o o ° 0 

aaaaaaa** A 


A A 


A A 


OO000000000 


A A 


X X 


X x 


X X 


o ° 


o 0 


-©OOOOOOOO 


0.0 


0.2 


0.4 


0.6 


0.8 


1 . 


4 > /*■ 


Figure 14. Boundar. -Layer Thickness (Cone, M ^ — 3.5, a — 2.25°) 
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Figure 15. Displacement Thickness (Cone, - 3.5, a - 2.25°) 
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Figure 17. Crosswise Velocity Profile at X = 0.582 ft (Cone, M — 3.5, a — 2.25 
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Figure 27. Sign ot m viscid Crossflow Velocity Component (Fuselage, Moo — 0-3, ct = 3°) 
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Figure 28. Inviscid Cross [low Velocity Component (u e ) as a Function of <j> 
(Fuselage, — 0.3, a — 3°) 
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Figure 35. Comparison of Three-Dimensional Velocity Profiles with Axisymmetric 
Analogue Results (Fuselage, M = 0.3, a = 3°, X = 1.3 m) 
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Figure 36. Skin Friction near Nose for Various Ax 

(Ellipsoid of Revolution, Fqo = 1 m/sec, a — 0°) 
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Figure 38. Streamline Coordinates Near Stagnation Point 
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Figure 39. Sketch :"-r Calculating the initial Velocity Profiles in the Streamline 


Coordinates Near Sharp Nose 




X-const line 



Figure 40. Sketch Showing the Velocity Components in Body-Oriented Coordinates 
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(b) Cross Section 


Figure 41. Spherical Polar Coordinates 
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